Cortical bone continuum damage mechanics constitutive model with stress triaxiality criterion to predict fracture initiation and pattern

A primary objective of finite element human body models (HBMs) is to predict response and injury risk in impact scenarios, including cortical bone fracture initiation, fracture pattern, and the potential to simulate post-fracture injury to underlying soft tissues. Current HBMs have been challenged to predict the onset of failure and bone fracture patterns owing to the use of simplified failure criteria. In the present study, a continuum damage mechanics (CDM) model, incorporating observed mechanical response (orthotropy, asymmetry, damage), was coupled to a novel phenomenological effective strain fracture criterion based on stress triaxiality and investigated to predict cortical bone response under different modes of loading. Three loading cases were assessed: a coupon level notched shear test, whole bone femur three-point bending, and whole bone femur axial torsion. The proposed material model and fracture criterion were able to predict both the fracture initiation and location, and the fracture pattern for whole bone and specimen level tests, within the variability of the reported experiments. There was a dependence of fracture threshold on finite element mesh size, where higher mesh density produced similar but more refined fracture patterns compared to coarser meshes. Importantly, the model was functional, accurate, and numerically stable even for relatively coarse mesh sizes used in contemporary HBMs. The proposed model and novel fracture criterion enable prediction of fracture initiation and resulting fracture pattern in cortical bone such that post-fracture response can be investigated in HBMs.

1 Introduction

The modeling of hard tissue response under load, fracture initiation, and post-fracture response for different modes of loading is critical for the prediction of Crash Induced Injuries in advanced finite element (FE) Human Body Models (HBMs) (Schmitt et al., 2019). Cortical bone is present in the diaphyses of long bones, as a thin shell in the epiphyses of long bones and surrounding the flat bones (Shore et al., 2013), where fractures involving cortical bone are important in injury biomechanics and may be associated with serious or greater injury levels (i.e., Abbreviated Injury Scale AIS 3+) (AAAM, 2015). Furthermore, predicting fracture pattern and post-fracture response of hard tissue is critical for the future possibility of assessing injury risk to underlying soft tissues. Current HBMs often model cortical bone using an isotropic, elastic-plastic material model incorporating a yield surface to predict the response of hard tissue (Yang, 2017), while failure is modeled using element erosion at a specified effective plastic strain (DeWit, & Cronin, 2012). Two benefits of this simplified approach are computational efficiency and numerical robustness, ensuring the calculation runs to completion with the coarse mesh sizes (e.g. 1–3 mm) used in contemporary HBMs. Such constitutive models can predict the onset of failure arising from uniaxial tension loading, since the material parameters and failure criterion are usually calibrated to this load case, but are limited in predicting failure initiation in other modes of loading (Khor et al., 2018), and have not been successful at predicting bone fracture pattern. Although detailed microstructural models have been proposed for cortical bone (Dapaah et al., 2020), such models with elements on the order of 10 microns in the region of interest, are too computationally expensive for whole body HBM simulations.

Current state-of-the-art HBMs require a constitutive model, which is functional for coarse element sizes (∼1–3 mm) used in HBMs, relative to contemporary micro-models. Previous research (Khor et al., 2018) has demonstrated the need for asymmetric tension-compression behavior and orthotropy. Cortical bone also exhibits damage or post-yield behavior and complex fracture patterns depending on the mode of loading (Figure 1). In the present study a continuum damage mechanics (CDM) material model with material direction-dependent damage was integrated with a proposed stress triaxiality-based fracture criterion to predict the onset of failure and resulting bone fracture pattern. The fracture criterion parameters were fit to experimental data. The resulting predicted fracture patterns were assessed using three experimental load cases.

FIGURE 1

FIGURE 1. Cortical bone tension response (A), microcrack damage (B), and resulting transverse fracture (C); shear response (D), microcrack damage (E), and resulting fracture (F).

1.1 Cortical bone structure, mechanical properties and constitutive models

Approximately 80% of the human skeletal bone mass is cortical bone (Cowin, 2001) comprising inorganic hydroxyapatite, organic Type I collagen phases, other collagen types, non-collagenous proteins, and water (Keaveny et al., 2003). The stiffness of bone has been associated with the inorganic phase (hydroxyapatite) (Saito, & Marumo, 2015); while the organic content (primarily Type I collagen) is associated with the tensile properties and post-yield deformation of the bone. The structure of cortical bone is hierarchical in nature, beginning with the inorganic and organic phases (∼0.1 micron), increasing to lamellae (∼1–10 micron) and cylindrically-shaped osteons (200 µm diameter, 1–3 mm length). Within the long bones of the body, the osteons are oriented along the diaphyseal or long axis of the bone (Idkaidek & Jasiuk, 2017). Cortical bone has a higher stiffness and strength along the longitudinal or osteon direction than the transverse direction, with orthotropic shear properties (Tang et al., 2015). In addition, the properties of cortical bone may vary with age, sex and lifestyle (Bruno et al., 2014), and the measured mechanical properties may vary with hydration, embalming (Sanborn et al., 2016; Öhman et al., 2008; Stefan et al., 2010) and specimen manufacturing method (Cowin, 2001). A previously reported set of mechanical properties was presented by Khor et al. (2018) (Table 1).

TABLE 1

TABLE 1. Cortical bone material property summary (Khor et al., 2018).

Under uniaxial tension loading, cortical bone material response is characterized in three phases (Figure 1A): an initial elastic response; continuum damage mechanics response, where damage in the material accumulates as microcracks oriented perpendicular to the applied load direction or principal stress (Figure 1B); and the fracture mechanics region where damage localizes leading to the initiation and propagation of a fracture (Gupta & Zioupos 2008). The microcracks occur perpendicular to the applied stress, providing energy absorption prior to failure, where the resulting crack paths typically follow the osteon cement lines (Zhai et al., 2019a). Experimental studies have reported sublaminar microcracking in bone but were focused on compressive loading (Ebacher et al., 2012), whereas the fracture of cortical bone has been reported to initiate in regions subjected to tensile loading (Zhai et al., 2019b; Wolfram & Schwiedrzik, 2016). Tension loading parallel to the osteons at the coupon level generally leads to transverse fracture of the specimen, perpendicular to the applied load (Figure 1C). The strain dependence of cortical bone failure has been reported (Nalla et al., 2003), owing to diffuse microcracking damage that may be pressure-dependent.

Several shear tests have been applied to test bone including the rail shear test, the torsion tube, and cross-beam specimen, with the limitation that some methods inadvertently generate tension in the test specimen (Cowin, 2001). Damage reported from shear loading (Figure 1D) occurs as microcracks in orthogonal directions (Figure 1E); however, failure is not reported to occur in the zone of maximum shear within the notched area or ligament of the test specimen. The Iosipescu test involves shear loading using rigid clamping fixtures to induce uniform shear in the ligament of the test specimen (Tang et al., 2015). Shear loading using the notched shear (Iosipescu) geometry results in a non-intuitive fracture initiation away from the notch and propagating at an angle into the sample (Figure 1F) (Tang et al., 2015), demonstrating the complex behavior of cortical bone.

Cortical bone constitutive models have been investigated at the micro-scale (Idkaidek & Jasiuk, 2017; Abdel-Wahab et al., 2012; Ascenzi et al., 2013; Demirtas et al., 2016; Ural & Vashishth, 2006; Ural et al., 2011; Feerick et al., 2013), at the macro-scale at coupon level (Garcia et al., 2009; Li et al., 2013), and at the scale of whole bones (San Antonio et al., 2012; Ariza et al., 2015; Asgharpour et al., 2014; Keyak et al., 1997; Iwamoto et al., 2005; Schileo et al., 2008; Zysset et al., 2013). Contemporary HBMs often use linear isotropic and symmetric metals plasticity constitutive models (DeWit, & Cronin, 2012; Asgharpour et al., 2014; Untaroiu et al., 2013). A recent study (Khor et al., 2018) identified the importance of including asymmetry and orthotropy to predict the response and failure of whole bones. In general, orthotropic material properties were essential to predict the onset of failure for whole bone simulations in two primary modes of loading (bending and shear); however, the predicted fracture patterns were not in agreement with reported patterns when using plasticity-based metal constitutive models. Incorporation of tension-compression asymmetry generally improved the gross kinetics and kinematics of whole bone fracture simulations, but similar to contemporary isotropic constitutive models, failure was predicted to initiate at the point of loading in an area of high compressive stress, which differed from the reported experiments where failure is reported in the tension region of the bone.

CDM models account for damage in a material by reducing the material stiffness to represent microcrack development. The Matzenmiller, Lubliner and Taylor (MLT) post-failure CDM (Matzenmiller et al., 1995; Gower et al., 2007) incorporates orthotropy, asymmetry and material damage. Within the model, damage (ω_ij) is represented as a scalar value from 0.0 to 1.0 for each material direction, with increasing damage corresponding to a reduction in material stiffness (E_ii or G_ij). Using Voight notation, the rate of damage accumulation,

$g=\left[\begin{array}{c}\begin{array}{c}{e}^{-1}\left(1-{\omega }_{11,\mathrm{t}-1}\right)\left(\frac{{\dot{\epsilon }}_{11}}{{\epsilon }_{11,\mathrm{f}}}\right){\left(\frac{{\epsilon }_{11}}{{\epsilon }_{11,\mathrm{f}}}\right)}^{{\mathrm{m}}_{11}-1}\\ e^{-1}\left(1-{\omega }_{22,\mathrm{t}-1}\right)\left(\frac{{\dot{\epsilon }}_{22}}{{\epsilon }_{22,\mathrm{f}}}\right){\left(\frac{{\epsilon }_{22}}{{\epsilon }_{22,\mathrm{f}}}\right)}^{{\mathrm{m}}_{22}-1}\\ e^{-1}\left(1-{\omega }_{33,\mathrm{t}-1}\right)\left(\frac{{\dot{\epsilon }}_{33}}{{\epsilon }_{33,\mathrm{f}}}\right){\left(\frac{{\epsilon }_{33}}{{\epsilon }_{33,\mathrm{f}}}\right)}^{{\mathrm{m}}_{33}-1}\end{array}\\ \begin{array}{c}{e}^{-1}\left(1-{\omega }_{12,\mathrm{t}-1}\right)\left(\frac{{\dot{\epsilon }}_{12}}{{\epsilon }_{12,\mathrm{f}}}\right){\left(\frac{{\epsilon }_{12}}{{\epsilon }_{12,\mathrm{f}}}\right)}^{{\mathrm{m}}_{12}-1}\\ e^{-1}\left(1-{\omega }_{23,\mathrm{t}-1}\right)\left(\frac{{\dot{\epsilon }}_{23}}{{\epsilon }_{23,\mathrm{f}}}\right){\left(\frac{{\epsilon }_{23}}{{\epsilon }_{23,\mathrm{f}}}\right)}^{{\mathrm{m}}_{23}-1}\\ e^{-1}\left(1-{\omega }_{31,\mathrm{t}-1}\right)\left(\frac{{\dot{\epsilon }}_{31}}{{\epsilon }_{31,\mathrm{f}}}\right){\left(\frac{{\epsilon }_{31}}{{\epsilon }_{31,\mathrm{f}}}\right)}^{{\mathrm{m}}_{31}-1}\end{array}\end{array}\right],(1)$

is defined by the damage exponent (m_ij) (Figure 2), which is fit to coupon-level material test data. Within the damage calculation, the current strain (ε_ij), current strain rate (ε̇_ij), material failure strain (ε_ij,f) defined as the material strength divided by the modulus, and damage parameter for a given material direction are used to determine the onset of damage through a damage initiation criterion. The damage rate multiplied by the time increment (dt) determines the increase in material damage for a given time increment,

${\omega }_{\mathrm{t}}={\omega }_{\mathrm{t}-1}+{\mathrm{g}}_t\mathrm{d}\mathrm{t},(2)$

and is applied in the material compliance tensor,

$H=\left[\begin{array}{cccccc}\frac{1}{\left(1-{\omega }_{11}\right)E_{11}}& \frac{-{\nu }_{21}}{E_{22}}& \frac{-{\nu }_{31}}{E_{33}}& 0& 0& 0\\ \frac{-{\nu }_{12}}{E_{11}}& \frac{1}{\left(1-{\omega }_{22}\right)E_{22}}& \frac{-{\nu }_{32}}{E_{33}}& 0& 0& 0\\ \frac{-{\nu }_{13}}{E_{11}}& \frac{-{\nu }_{23}}{E_{22}}& \frac{1}{\left(1-{\omega }_{33}\right)E_{33}}& 0& 0& 0\\ 0& 0& 0& \frac{1}{\left(1-{\omega }_{12}\right)G_{12}}& 0& 0\\ 0& 0& 0& 0& \frac{1}{\left(1-{\omega }_{23}\right)G_{23}}& 0\\ 0& 0& 0& 0& 0& \frac{1}{\left(1-{\omega }_{31}\right)G_{31}}\end{array}\right](3)$

FIGURE 2

FIGURE 2. CDM orthotropic and asymmetric model showing damage exponent (m) effect in uniaxial tension, and calibration for cortical bone tensile data (m = 0.726)

Damage coupling between directions is included, and is visible in the inverted compliance (stiffness) form of the matrix. Importantly, the CDM approach can distinguish between tension and compression and therefore include material asymmetry. The stresses corresponding to the damaged material state are calculated considering the material damage,

$\sigma =H^{-1}\epsilon .(4)$

Fracture initiation and propagation in bone has been investigated using cohesive-zone element models (Ural & Vashishth, 2006; Ural et al., 2011) with the limitation that the crack-extension path must be predefined (Li et al., 2013). Although the extended finite element method (X-FEM) allows for crack propagation to be modelled without the need to predefine the crack path (Li et al., 2013), current cortical bone models utilizing the X-FEM method are only two-dimensional and are often applied at the micro-scale level (Budyn et al., 2008; Feerick et al., 2013; Abdel-Wahab et al., 2012; Idkaidek & Jasiuk, 2017) prohibiting general use in HBMs.

It has been demonstrated that hydrostatic stress state plays a role in the failure of low ductility materials (Johnson & Holmquist, 1992; Morin et al., 2010). Recent developments in modeling fracture have identified the dependence of failure parameters on stress triaxiality (Butcher & Abedini 2019) for fracture of metals. Triaxiality (η) is defined as the ratio of hydrostatic stress to effective stress, such that a triaxiality of zero corresponds to pure shear loading and a value of 1/3 corresponds to uniaxial tension. The element deletion method (element erosion) with a strain-based failure criterion is still widely used to predict material fracture, and has been somewhat successful at predicting the onset of failure (DeWit, & Cronin, 2012; Schileo et al., 2008; Niebur et al., 2000). The element erosion method is often numerically stable, but has generally been unable to predict bone fracture patterns for various modes of loading.

1.2 Experimental testing of whole bones

Whole bone testing is considered to be the most relevant representation of bone response and failure (Shore et al., 2013) since there are no specimen machining artifacts induced. Several studies have been undertaken to measure the mechanical response of whole bones in the lower extremity (Cezayirlioglu et al., 1985; Kress et al., 1995; Strømsøe et al., 1995; Cristofolini et al., 1996; Funk et al., 2004). Khor et al. (2018) identified two whole bone femur load cases to assess constitutive models: a three-point bending load case, including both posterior-anterior and medial-lateral loading conditions (Funk et al., 2004; Kerrigan et al., 2003), and an axial torsion load case (Martens et al., 1980), using an apparatus designed by Burstein and Frankel (1971). The experimental results (Table 2) demonstrated expected fracture patterns for torsional and bending loading.

TABLE 2

TABLE 2. Whole bone femur experimental data.

Initiation of cortical bone fracture is often associated with tensile loading, and locally a defect or stress concentration on the bone surface determines the fracture onset location (Carter & Spengler, 1982). Following initiation, the fracture propagates approximately perpendicular to the maximum principal stress. In vivo fractures often exhibit complex fracture patterns due to complex loading (Carter & Spengler, 1982) and higher loading rates tend to increase comminution or fragmentation of the fracture (Kress et al., 1995; Turner, 2006). However, single mode loading fracture patterns are well-established. The strength of cortical bone is lowest in tension and shear, compared to compression (Turner, 2006) and therefore cracks tend to propagate along tension or shear planes within the bone tissue. The fracture pattern associated with tensile loading is a transverse fracture oriented perpendicular to the longitudinal or osteon direction. In bending, failure initiates on the tensile side of the bone due to the lower strength in tension compared to compression (Kress et al., 1995; Turner, 2006; Sharir et al., 2008). As the crack propagates transversely across the bone towards the compressive stress zone, the fracture often bifurcates at an angle of approximately 45° (Turner, 2006; Sharir et al., 2008). This is known as a butterfly or tension wedge fracture and is a typical fracture pattern for bending loading (Sharir et al., 2008). Torsion loading results in a spiral fracture pattern, explained by the maximum tension plane located at 45° from the shear plane, causing the crack to propagate along a helical plane of maximum tension (Turner, 2006).

2 Methods

The present study investigated cortical bone response using three models: a contemporary metals plasticity model (Khor et al., 2018), a CDM model using measured mechanical properties of cortical bone with fracture occurring based on damage accumulation, and the CDM model including a novel phenomenological fracture criterion. Bone fracture patterns were assessed using three test cases: notched shear (Iosipescu) specimen, whole femur three-point bending, and whole femur axial torsion. The models were solved in a commercial explicit FE program (LS-DYNA R9.3, LST, Livermore, CA) compiled with a custom code (user material model) for the CDM model and fracture criterion, described below.

2.1 Iosipescu (notched shear) specimen simulations

Models of the notched shear test specimen (Figure 3, blue) and loading apparatus (grey) were created in a commercial preprocessor software (Hypermesh, Altair) using the specimen dimensions and boundary conditions described in the experiments (Tang et al., 2015). The loading apparatus was modeled with an elastic material with the properties of steel, and loading was applied at a rate of 0.01 m/s using a displacement boundary condition as in the experiments. The notched shear simulation was used initially to investigate failure criteria and to develop the proposed failure criterion. The model response was verified by comparing to the experimental force-displacement response.

FIGURE 3

FIGURE 3. Finite element models of load cases and boundary conditions.

2.2 Whole femur finite element model

The femur finite element model was extracted from a current average stature (50th percentile) male human body model (M50 Version 6.0, Global Human Body Models Consortium). In a previous study, this femur model was shown to fall within one standard deviation in length and cross-sectional area reported in the experiments (Khor et al., 2018). Owing to the continuum nature of the constitutive model in this study, the osteons were not explicitly modeled; however, the effect of the osteon orientation and connectivity is known to be associated with fracture patterns and was included in the model by assigning the nodes of each element to define element-level material directions (Supplementary Figure S1). The primary material direction was the osteon direction, and the two orthogonal directions were radial (through thickness) and circumferential, as required for the orthotropic CDM model.

2.2.1 Three-Point bending whole femur load case

The three-point bending boundary conditions were applied to the femur at the epiphyses (Figure 3). The proximal femur translation in the model was fixed in the inferior-superior, anterior-posterior, and medial-lateral directions (x, y and z-directions, respectively, in the model). The distal femur translation was fixed in the anterior-posterior, and medial-lateral directions (y and z-directions, respectively, in the model). Both epiphyses were free to rotate about the direction perpendicular to the loading, as described in the experimental studies. The long axis of the bone was aligned with the global x-axis of the model and loading was applied by a steel semi-cylinder near the midpoint along the diaphysis, corresponding to the loading point in the experimental tests. The semi-cylinder was moved (1.2 m/s) using a prescribed displacement condition and load was monitored through the contact force. The contact algorithm between the semi-cylinder and the bone accounted for element erosion, and reformulated the contact surface after each eroded element was removed to represent the loading and interaction with the bone as fracture progressed. The model response was validated by comparing to the experimental force-displacement response.

2.2.2 Axial Torsion whole femur load case

The femur was aligned with the global y-direction for the axial torsion load case (Figure 3) with the distal epiphysis of the bone fixed in rotation and translation, and the proximal epiphysis loaded using a prescribed rotation displacement (0.00873 rad/ms or 500^°/s). The model response was validated by comparing to the experimental torque-rotation response.

2.3 Orthotropic continuum damage mechanics model and triaxiality-based failure criterion

Within the current study, three material models were analyzed. The reference case was an isotropic plasticity model as implemented in contemporary HBMs, with failure predicted by element erosion at an effective plastic strain of 1.8% (Khor et al., 2018). The second case was the MLT CDM, implemented with element failure occurring at a damage value of 1.0. The MLT implementation enabled for the assessment of material orthotropy, tension-compression asymmetry and damage. Thirdly, a novel stress triaxiality failure criterion was integrated with the MLT CDM model, described as the Cortical Bone Fracture and Continuum Damage Mechanics Model (CFraC).

2.3.1 MLT continuum damage mechanics model

The MLT CDM model was implemented as a user-defined material model as reported in the literature (Gower et al., 2007), in a commercial explicit FE program. The model enabled simulation of the elastic and the CDM portions of the tensile response (Figure 1) using published material properties. The material model parameters were determined as follows:

1) Elastic moduli, Poisson’s ratios and strength values were taken from the literature (Table 1). Although cortical bone properties do exhibit variability, the average elastic mechanical properties and strength properties are generally agreed upon in the literature (Khor et al., 2018).

2) The tensile damage exponent in the osteon direction (1-direction) (m₁₁) and transverse directions (m₂₂, m₃₃) were determined using single element simulations of uniaxial tensile test data in the respective directions, with equivalent strain energy density in the model and experiment (i.e. the same area under the stress-strain curve) (Figure 2; Supplementary Figure S2). In general, the shape of the MLT CDM model differed from the reported experimental data, in that the experimental data demonstrated a more abrupt change from elastic to CDM regions, with a shallow slope in the CDM region. This difference is a limitation of the MLT model formulation and should be investigated for future CDM approaches.

3) It was assumed that damage and failure at the element level did not occur directly in compression or shear (i.e., ω₁₂, ω₂₃ and ω₃₁ were zero in Eq. 3). However, failure of an element could occur in these modes of loading due to the implemented damage coupling with other material directions.

2.3.2 Stress triaxiality-based failure criterion (CFraC)

The observation that hydrostatic stress plays a role in the failure of low ductility materials, and the dependency of failure on stress triaxiality led to the investigation of a novel criterion to represent the fracture mechanics (failure) region of cortical bone (Figure 1A) in the present study. It was hypothesized that damage accumulation and localization may be a function of stress triaxiality and effective strain. This hypothesis, which still requires experimental verification, was pursued in the current study. The proposed fracture criterion was coupled to the MLT CDM model to create the Cortical Bone Fracture and Continuum Damage Mechanics Model (CFraC) and investigated using three test cases (Iosipescu test, three-point bending, axial torsion). A fourth case (uniaxial tension) was introduced with a known tensile hydrostatic stress (1/3 of the material strength in a given direction) and stress triaxiality (1/3). The effective strain versus triaxiality curve was determined by simulating each of the three load cases (Iosipescu test, three-point bending, axial torsion) without failure. At the force or torque corresponding to failure in the experiments, the effective strain and corresponding triaxiality were determined for each load case. Element failure (erosion) was based on the effective strain versus triaxiality curve constructed using each of the load cases. The critical hydrostatic stress values were determined from the finite element femur model so that the presented values were relevant to the finite element mesh sizes used in contemporary HBMs. However, finite element models and simulation of failure processes are known to have a dependence on finite element mesh size. In the present study, these effects were investigated using two refined models with mesh densities increased by a factor of 2x and 4x. The results were compared between different mesh sizes to assess the effect of mesh refinement on the predicted mechanical response and failure.

To verify the material model implementation in the FE code, single element simulations were undertaken for the osteon and transverse directions in tension, compression and shear. A mesh size of 3 mm was used for these simulations. In each case, the element was fixed in the loading direction with one-eighth symmetry conditions applied, and displacement boundary conditions were applied to the opposite element face. The element was not constrained opposite the symmetry planes, allowing for deformation due to Poisson’s ratio effects. Implementing the orthotropic material properties (Table 1) in the MLT model, single element simulations were undertaken in longitudinal (osteon direction) tension/compression, orthogonal (circumferential and radial directions) tension/compression, and in-plane shear direction. All models provided the expected stress-strain response, which verified the model implementation and input data.

3 Results and discussion

3.1 Failure criterion: Effective strain versus stress triaxiality

The fracture initiation threshold curve, describing hydrostatic stress as a function of triaxiality (η), was plotted for the three load cases (data points, Figure 4 and Table 3), along with the reported uniaxial tensile failure value at a triaxiality of 1/3. The data were then fit to a Modified Mohr-Coulomb model (Kõrgesaar et al., 2018), which was implemented in the computational model to define the effective strain at failure (ε̅_f) versus triaxiality:

${\overline{\epsilon }}_f={\left\{\frac{K_m}{C_2}{f}_3\left[\sqrt{\frac{1+C_1^2}{3}}{f}_1+C_1\left(\eta +\frac{f_2}{3}\right)\right]\right\}}^{-1/n_m}(\mathrm{5a})$

$f_1=\cos \left\{\frac{1}{3}\arcsin \left[-\frac{27}{2}\eta \left({\eta }^2-\frac{1}{3}\right)\right]\right\},(\mathrm{5b})$

$f_2=\sin \left\{\frac{1}{3}\arcsin \left[-\frac{27}{2}\eta \left({\eta }^2-\frac{1}{3}\right)\right]\right\},(\mathrm{5c})$

$f_3=C_3+\frac{\sqrt{3}}{2-\sqrt{3}}\left(1-C_3\right)+\left(\frac{1}{f_1}-1\right).(\mathrm{5d})$

FIGURE 4

FIGURE 4. Effective strain versus stress triaxiality (η) failure criterion.

TABLE 3

TABLE 3. Modified Mohr-Coulomb fracture model parameters.

Although the proposed failure locus is empirical at present, it highlights potential need for additional experimental information considering varying triaxiality, as was achieved in the Iosipescu samples (Tang et al., 2015) with high triaxiality at the failure location adjacent to the notch.

3.2 Iosipescu notched shear simulation

The notched shear test simulation demonstrated a monotonically increasing force versus displacement response up to the initiation of failure (Figure 5). The isotropic plasticity model significantly over predicted the failure force in shear due to the yield surface assumption in the model, highlighting a significant limitation of plasticity approaches to modeling cortical bone under shear loading.

FIGURE 5

FIGURE 5. Force-displacement response for notched shear (Iosipescu) specimen (experimental response from Tang et al., (2015) (A); and predicted fracture via eroded elements (B), fracture initiation location indicated by arrow.

Although the MLT CDM model predicted the failure force (0.5 kN) in reasonable agreement with the experiments, owing to input of the experimentally reported shear strength, failure was predicted in the central, shear gauge area between the two notches of the specimen, which disagreed with the experimental data (Supplementary Figure S3).

The CFraC model predicted a failure force (−0.41 kN), in agreement with the average reported experimental value (−0.42 kN). The simulation exhibited a similar stiffness to the experiment, but did not incorporate a toe region, attributed to compliance in the test fixture and the challenge of experimentally measuring very small displacements. The Iosipescu test specimen demonstrated high triaxiality (η = 0.54) at the fracture initiation location, approaching equibiaxial tension (η = 0.666). Fracture was predicted to initiate adjacent to the root of the sample notch and propagated at an approximate angle of approximately 40° to the horizontal plane, in agreement with the experimental findings.

3.3 Whole femur three-point bending simulation

For the whole femur three-point bending load case, the model force-displacement response monotonically increased up to the initiation of failure (Figure 6). The isotropic plasticity model failure force was below the test average and range owing to failure initiation at the load point early in the load history. Progressive crushing of the bone at the load point occurred providing an apparent energy absorption, but the fracture pattern did not agree with those reported in the literature.

FIGURE 6

FIGURE 6. Force-displacement response for whole bone posterior-anterior bending test (experimental response from Funk et al., (2004) (A); and predicted tension-wedge fracture via eroded elements (B), fracture initiation location indicated by arrow.

The MLT model predicted a failure force above to the experimental average. Although the fracture initiated on the tension side of the bone, the fracture bifurcated to follow the neutral axis of the bone (Supplementary Figure S3) due to the lower shear strength of the bone, and was not representative of reported fracture patterns.

The maximum force of the CFraC model in posterior-anterior bending (4.15 kN) was within the range of the experimental data (3.6–4.9 kN). The displacement at failure (14.6 mm) agreed with the experimental average (16.7 mm, standard deviation of 3.38 mm). In medial-lateral bending, the maximum force was 3.5 kN at a displacement of 15.8 mm. Additional simulations demonstrated that achieving a higher force value required an increase in the tensile material strength (Table 1) or decreasing the damage exponent (m11). The tensile strength of the bones tested in the experiments conducted by Funk et al. (2004) may have differed from properties reported by Reilly et al. (1974) that were used in the constitutive model. The geometry of the femur used in the present study (length and cross-sectional area) was similar to the average reported in the experiments (Khor et al., 2018). However, variations within the population are expected to result in a range of physical dimensions and material properties, which may explain some of the variability in the experimental test results (Table 2). The nonlinear response of the MLT and CFraC models was due to the accumulation of damage and corresponding reduction in stiffness of the cortical bone material. The CFraC model predicted failure initiation on the tension surface of the bone and fracture propagated across the bone to the compression side, as reported in the literature.

3.4 Whole femur axial torsion simulation

When the femur was loaded in axial torsion, the moment rotation response was monotonically increasing and relatively linear up to the point of failure. The isotropic plasticity model predicted a very high torque at failure for the whole femur axial torsion load case (Figure 7), attributed to the yield surface definition in the material model. The MLT model also predicted a torque and axial rotation exceeding the reported experimental range with a linear fracture pattern along the length of the bone (Supplementary Figure S3), rather than a spiral fracture, owing to the low shear strength along the length of the bone.

FIGURE 7

FIGURE 7. Torque-rotation response for whole bone axial torsion test (experimental response from Martens et al., (1980) (A); and predicted spiral fracture via eroded elements (B), fracture initiation location indicated by orange arrow.

For the CFraC model, failure initiated at a torque of 182 Nm in the proximal end of the diaphysis, and propagated in a spiral fracture towards the distal end of the bone. The predicted rotation to failure (19°) was within the experimental bounds (9.4–30.7°) and comparable to the reported experimental average (20°). The axial rotation for a given axial torque was determined by the torque, polar second moment of area, femur effective length, and shear modulus. Thus, increased length, reduced cross-sectional area or reduced stiffness could affect the rotation but were not investigated in the current study. While triaxiality is zero for pure shear loading, the triaxiality in the diaphysis under axial torsion ranged from zero up to ∼0.1. The varying η value was associated with changes in geometry of the bone. It should be noted that the failure location was within the diaphysis, and away from the epiphyses and applied boundary condition effects. This small but non-zero η value highlights the potential importance of the local bone geometry on response and failure prediction.

3.5 CFraC finite element model mesh refinement

A mesh refinement study was undertaken for the CFraC constitutive and fracture model using the three load cases. Finite element models are known to have a dependence on finite element size, attributed in part to the modeling of stress or strain gradients, while modeling of failure may be a non-convergent phenomenon. The finite element meshes were refined by splitting all solid elements (1 element split to 8 elements) and all shell elements (1 element split to 4 elements). This procedure was carried out twice to provide two refined meshes per model, in addition to the baseline mesh. All analyses were run to the same termination time, using the same boundary conditions and material properties. In general, refining the mesh resulted in a slightly lower predicted failure force or moment and failure displacement or rotation (Supplementary Figure S4), due to the improved resolution of the high hydrostatic tension gradient in the failure zone. Each refined model demonstrated the same location of fracture initiation and a more refined fracture pattern. Since the objective of this study was to provide a cortical bone material model for use in contemporary HBM, the critical effective strain was determined from the coarse mesh (production mesh for the GHBMC M50 model).

3.6 General discussion

The proposed failure initiation criterion was fit to the experimental data using the simulations so it should be expected that the models would predict the onset of failure in agreement with the experiments. In practice, mechanical tests are undertaken using specimen geometries that achieve varying levels of triaxiality using, for example, tensile test specimens with various notch geometries. The Iosipescu sample presents an example of this type of test. Although the goal of this shear test was to achieve pure shear loading (η = 0) at the failure location, this case demonstrated the highest triaxiality of all cases considered. In general, material failure data as a function of triaxiality is not currently available for cortical bone; however, the current model suggest this experimental data would be useful to further refine the predictive capabilities of cortical bone models.

One aspect not considered in the current model are deformation rate effects. Studies have identified deformation rate effects in compression loading (McElhaney, 1966) comprising increased modulus, increased strength and decreased strain to failure with increasing loading rates. However, studies are generally not in agreement regarding deformation rate effects in tension loading, where some researchers have measured increasing (Melnis & Knets, 1982) and reduced (Hansen et al., 2008) stiffness and strength with increasing loading rates. Given that tension appears to be an important factor in accumulating damage, and experimental studies to date present inconsistent information, deformation rate effects were not included in the current study. Further mechanical testing is needed to better quantify the effect of deformation rate on the mechanical properties of bone. The present study investigated fractures as reported in the literature that were generated from dynamic loading relevant to automotive crash or sports injuries. Additional three-point bending and torsion testing across a range of loading rates would provide some insight into the evolution of the facture pattern with loading rate. This potential evolution could then be considered in the model by implementing a strain-rate dependent fracture criterion. It is important to note that very high rate loading such as that encountered in blast exposure or ballistic impact may result in less well-defined fracture patterns, and more comminuted fracture. For example, exposure to antipersonnel landmines often results in comminuted fractures of foot and long bones (Cronin et al., 2011) while ballistic impacts on bone may result in drill-hole type fractures in the ribs and sternum (Nguyen et al., 2022).

The proposed model was symmetric in shear (in-plane direction) while the bone structure is asymmetric in the two in-plane shear directions owing to the orientation of the osteons and cement lines. Therefore, only shear across the osteons or parallel to the osteons can be represented. Within the current model, shear across the osteons was incorporated into the model, and provided good predictions of fracture strength and pattern for the axial torsion and Iosipescu cases. Material data in shear is somewhat limited, owing to the challenges in experimentally achieving pure shear loading, and the complex failure modes exhibited by cortical bone. The study by Tang et al. (2015) demonstrated that shear loading of cortical bone transverse to the osteon direction leads to complex failure behavior, although the measured force-displacement response is believed to represent the material shear behavior up to the point of fracture initiation. The Iosipescu sample should have zero triaxiality through the intended gauge section, but actually demonstrates very high triaxiality at the failure location owing to the presence of the notch. Current work in fracture of materials uses various notched samples in tension and shear to measure material failure for different levels of η. Such testing should be undertaken for bone in the future. The present study provides an estimate of the failure curve by interpreting existing experiments, but was limited in terms of material data available for a range of η values.

In the models of the Iosipescu tests, the displacement to failure was under predicted compared the experiments. This may be related to compliance in the experimental test apparatus, and variations between the moduli used in the simulations compared to the tests.

The constitutive and fracture model was computationally stable in all cases, allowing for the prediction of loading, fracture, and post-fracture response of the bone. Importantly, the CFraC model was able to predict the ultimate load, initiation of fracture and fracture pattern for the common modes of loading. The computation time for the whole femur subjected to axial torsion was 1.59 times greater for the CFraC model (Supplementary Figure S5) compared to the isotropic plasticity model, when simulated on one compute core of a Symmetric Multiprocessing system (i9 9960x CPU at 3.1 GHz). This ratio decreased to 1.02 when using 16 cores. In general, the CFraC model was more computationally expensive than the isotropic plasticity, which is anticipated owing to the larger number of calculations required for the CFraC model.

4 Study conclusions and limitations

The isotropic material model used in contemporary HBMs was able to predict the failure force associated with tension-based failures (e.g. three-point bending) but significantly over predicted the failure force or torque for shear loading due to the yield surface assumption embedded in metals plasticity models. An orthotropic MLT CDM material model including tension-compression asymmetry provided improved kinetic predictions for the load cases considered but was not able to predict fracture patterns in agreement with reported observations.

The proposed fracture initiation criterion (CFraC) was phenomenological, based on available experimental data and calibrated to the whole bone tests with respect to the failure initiation threshold; however, the resulting fracture patterns were not calibrated in any way. Importantly, the proposed failure criterion enabled the prediction of cortical bone fracture patterns in agreement with experimentally observed data. Refining the finite element mesh resulted in reduced force or moment and displacement to failure, as expected; and led to more refined fracture patterns that were in good agreement with reported bone fracture patterns for different modes of loading.

The study employed a subject-specific femur bone model, with length and cross-sectional area similar to the average values presented in experimental studies, and therefore was expected to be representative for the purposes of the current study. Variability in material properties and femur geometry were not investigated in the current study but may explain some of the differences between the model predictions and the average experimental results.

The CFraC constitutive and fracture model was computationally stable in all three load cases, allowing for the prediction of loading, fracture, and post-fracture response of the bone. Importantly, the model was able to predict the ultimate load, initiation of fracture and fracture pattern for the common modes of loading.

Data availability statement

The novel material model developed in this study will be available to the research community as a shared library that can be linked into the LS-DYNA finite element code. Please contact the corresponding author for information on accessing the material model.

Author contributions

DC implemented the material model in the finite element code and simulated the various boundary conditions. BW, FK, DG and SM developed models, assisted in data analysis, and generating figures for the manuscript. All authors contributed to this manuscript and approved the submitted version.

Funding

Fund for the study was provided by Honda Development and Manufacturing of America.

Acknowledgments

The authors would like to express their appreciation to Honda Development and Manufacturing of America for financial support of this research, the Global Human Body Models Consortium for use of the human model, and the Digital Research Alliance of Canada for providing the necessary computing resources.

Conflict of interest

Author SM was employed by the company Honda Development and Manufacturing of America.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2022.1022506/full#supplementary-material

References

AAAM (Association for the Advancement of Automotive Medicine) (2015). Abbreviated injury scale. Available at: https://www.aaam.org/abbreviated-injury-scale-ais/.

Google Scholar

Abdel-Wahab, A. A., Maligno, A. R., and Silberschmidt, V. V. (2012). Micro-scale modelling of bovine cortical bone fracture: Analysis of crack propagation and microstructure using X-FEM. Comput. Mater. Sci. 52 (1), 128–135. doi:10.1016/j.commatsci.2011.01.021

CrossRef Full Text | Google Scholar

Ariza, O., Gilchrist, S., Widmer, R. P., Guy, P., Ferguson, S. J., Cripton, P. A., et al. (2015). Comparison of explicit finite element and mechanical simulation of the proximal femur during dynamic drop-tower testing. J. biomechanics 48 (2), 224–232. doi:10.1016/j.jbiomech.2014.11.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Ascenzi, M. G., Kawas, N. P., Lutz, A., Kardas, D., Nackenhorst, U., and Keyak, J. H. (2013). Individual-specific multi-scale finite element simulation of cortical bone of human proximal femur. J. Comput. Phys. 244, 298–311. doi:10.1016/j.jcp.2012.05.027

CrossRef Full Text | Google Scholar

Asgharpour, Z., Zioupos, P., Graw, M., and Peldschus, S. (2014). Development of a strain rate dependent material model of human cortical bone for computer-aided reconstruction of injury mechanisms. Forensic Sci. Int. 236, 109–116. doi:10.1016/j.forsciint.2013.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Ashman, R. B., Cowin, S. C., Van Buskirk, W. C., and Rice, J. C. (1984). A continuous wave technique for the measurement of the elastic properties of cortical bone. J. biomechanics 17 (5), 349–361. doi:10.1016/0021-9290(84)90029-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruno, A. G., Broe, K. E., Zhang, X., Samelson, E. J., Meng, C. A., Manoharan, R., et al. (2014). Vertebral size, bone density, and strength in men and women matched for age and areal spine BMD. J. Bone Min. Res. 29 (3), 562–569. doi:10.1002/jbmr.2067

PubMed Abstract | CrossRef Full Text | Google Scholar

Budyn, E., Hoc, T., and Jonvaux, J. (2008). Fracture strength assessment and aging signs detection in human cortical bone using an X-FEM multiple scale approach. Comput. Mech. 42 (4), 579–591. doi:10.1007/s00466-008-0283-1

CrossRef Full Text | Google Scholar

Burstein, A. H., and Frankel, V. H. (1971). A standard test for laboratory animal bone. J. biomechanics 4 (2), 155–158. doi:10.1016/0021-9290(71)90026-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Butcher, C., and Abedini, A. (2019). On phenomenological failure loci of metals under constant stress states of combined tension and shear: Issues of coaxiality and non-uniqueness. Metals 9 (10), 1052. doi:10.3390/met9101052

CrossRef Full Text | Google Scholar

Carter, D. R., and Spengler, D. M. (1982). “Biomechanics of fracture,” in Bone in clinical orthopedics. Editor G. Sumner-Smith (Philadelphia: WB Saunders), 305–334.

Google Scholar

Cezayirlioglu, H., Bahniuk, E., Davy, D. T., and Heiple, K. G. (1985). Anisotropic yield behavior of bone under combined axial force and torque. J. biomechanics 18 (1), 61–69. doi:10.1016/0021-9290(85)90045-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Cowin, S. C. (2001). Bone mechanics handbook. 2nd Ed. Boca Raton: CRC Press. doi:10.1201/b14263

CrossRef Full Text | Google Scholar

Cristofolini, L., Viceconti, M., Cappello, A., and Toni, A. (1996). Mechanical validation of whole bone composite femur models. J. biomechanics 29 (4), 525–535. doi:10.1016/0021-9290(95)00084-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Cronin, D. S., Williams, K. V., and Salisbury, C. (2011). Physical surrogate leg to evaluate blast mine injury. Mil. Med. 176 (12), 1408–1416. PMID: 22338357. doi:10.7205/milmed-d-11-00044

PubMed Abstract | CrossRef Full Text | Google Scholar

Dapaah, D., Badaoui, R., Bahmani, A., Montesano, J., and Willett, T. (2020). Modelling the micro-damage process zone during cortical bone fracture. Eng. Fract. Mech. 224, 106811. doi:10.1016/j.engfracmech.2019.106811

CrossRef Full Text | Google Scholar

Demirtas, A., Curran, E., and Ural, A. (2016). Assessment of the effect of reduced compositional heterogeneity on fracture resistance of human cortical bone using finite element modeling. Bone 91, 92–101. doi:10.1016/j.bone.2016.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

DeWit, J. A., and Cronin, D. S. (2012). Cervical spine segment finite element model for traumatic injury prediction. J. Mech. Behav. Biomed. Mater. 10, 138–150. doi:10.1016/j.jmbbm.2012.02.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebacher, V., Guy, P., Oxland, T. R., and Wang, R. (2012). Sub-lamellar microcracking and roles of canaliculi in human cortical bone. Acta biomater. 8 (3), 1093–1100. doi:10.1016/j.actbio.2011.11.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Feerick, E. M., Liu, X. C., and McGarry, P. (2013). Anisotropic mode-dependent damage of cortical bone using the extended finite element method (XFEM). J. Mech. Behav. Biomed. Mater. 20, 77–89. doi:10.1016/j.jmbbm.2012.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Funk, J. R., Kerrigan, J. R., and Crandall, J. R. (2004). “Dynamic bending tolerance and elastic-plastic material properties of the human femur,” in Annual proceedings/association for the advancement of automotive medicine, 48, 215.

PubMed Abstract | Google Scholar

Garcia, D., Zysset, P. K., Charlebois, M., and Curnier, A. (2009). A three-dimensional elastic plastic damage constitutive law for bone tissue. Biomech. Model. Mechanobiol. 8 (2), 149–165. doi:10.1007/s10237-008-0125-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Gower, H. L., Cronin, D. S., and Plumtree, A. (2007). Ballistic impact response of laminated composite panels. Int. J. Impact Eng. 35 (9), 1000–1008. doi:10.1016/j.ijimpeng.2007.07.007

CrossRef Full Text | Google Scholar

Gupta, H. S., and Zioupos, P. (2008). Fracture of bone tissue: The ‘hows’ and the ‘whys. Med. Eng. Phys. 30 (10), 1209–1226. doi:10.1016/j.medengphy.2008.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, U., Zioupos, P., Simpson, R., Currey, J. D., and Hynd, D. (2008). The effect of strain rate on the mechanical properties of human cortical bone. J. Biomech. Eng. 130 (1), 011011. doi:10.1115/1.2838032

PubMed Abstract | CrossRef Full Text | Google Scholar

Idkaidek, A., and Jasiuk, I. (2017). Cortical bone fracture analysis using XFEM–case study. Int. J. Numer. Method. Biomed. Eng. 33 (4), e2809. doi:10.1002/cnm.2809

PubMed Abstract | CrossRef Full Text | Google Scholar

Iwamoto, M., Miki, K., and Tanaka, E. (2005). Ankle skeletal injury predictions using anisotropic inelastic constitutive model of cortical bone taking into account damage evolution. Stapp Car Crash J. 49, 133–156. doi:10.4271/2005-22-0007

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, G. R., and Holmquist, T. J. (1992). “A computational constitutive model for brittle materials subjected to large strains,” in Shock-wave and high strain-rate phenomena in materials. Editors M. A. Meyers, L. E. Murr, and K. P. Staudhammer (New York: Marcel Dekker), 1075–1081.

Google Scholar

Keaveny, T. M., Morgan, E. F., and Yeh, O. C. (2003). “Bone mechanics,” in Standard handbook of biomedical engineering and design. Editor M. Kurtz (New York: McGraw-Hill), 8.

Google Scholar

Kerrigan, J. R., Bhalla, K. S., Madeley, N. J., Funk, J. R., Bose, D., and Crandall, J. R. (2003). Experiments for establishing pedestrian-impact lower limb injury criteria. doi:10.4271/2003-01-0895

CrossRef Full Text | Google Scholar

Keyak, J. H., Rossi, S. A., Jones, K. A., and Skinner, H. B. (1997). Prediction of femoral fracture load using automated finite element modeling. J. biomechanics 31 (2), 125–133. doi:10.1016/S0021-9290(97)00123-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Khor, F., Cronin, D. S., Watson, B., Gierczycka, D., and Malcolm, S. (2018). Importance of asymmetry and anisotropy in predicting cortical bone response and fracture using human body model femur in three-point bending and axial rotation. J. Mech. Behav. Biomed. Mater. 87, 213–229. doi:10.1016/j.jmbbm.2018.07.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Kõrgesaar, M., Romanoff, J., Remes, H., and Palokangas, P. (2018). Experimental and numerical penetration response of laser-welded stiffened panels. Int. J. Impact Eng. 114, 78–92. doi:10.1016/j.ijimpeng.2017.12.014

CrossRef Full Text | Google Scholar

Kress, T. A., Porta, D. J., Snider, J. N., Fuller, P. M., Psihogios, J. P., Heck, W. L., et al. (1995). “Fracture patterns of human cadaver long bones,” in Proceedings of the International Research Council on the Biomechanics of Injury conference, 155–169.23

Google Scholar

Li, S., Abdel-Wahab, A., and Silberschmidt, V. V. (2013). Analysis of fracture processes in cortical bone tissue. Eng. Fract. Mech. 110, 448–458. doi:10.1016/j.engfracmech.2012.11.020

CrossRef Full Text | Google Scholar

Martens, M., Van Audekercke, R., De Meester, P., and Mulier, J. C. (1980). The mechanical characteristics of the long bones of the lower extremity in torsional loading. J. biomechanics 13 (8), 667–676. doi:10.1016/0021-9290(80)90353-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Matzenmiller, A. L. J. T. R., Lubliner, J., and Taylor, R. L. (1995). A constitutive model for anisotropic damage in fiber-composites. Mech. Mater. 20 (2), 125–152. doi:10.1016/0167-6636(94)00053-0

CrossRef Full Text | Google Scholar

McElhaney, J. H. (1966). Dynamic response of bone and muscle tissue. J. Appl. physiology 21 (4), 1231–1236. doi:10.1152/jappl.1966.21.4.1231

PubMed Abstract | CrossRef Full Text | Google Scholar

Melnis, A. E., and Knets, I. V. (1982). Effect of the rate of deformation on the mechanical properties of compact bone tissue. Mech. Compos. Mat. 18 (3), 358–363. doi:10.1007/BF00604319

CrossRef Full Text | Google Scholar

Morin, D., Haugou, G., Bennani, B., and Lauro, F. (2010). Identification of a new failure criterion for toughened epoxy adhesive. Eng. Fract. Mech. 77 (17), 3481–3500. doi:10.1016/j.engfracmech.2010.09.016

CrossRef Full Text | Google Scholar

Nalla, R., Kinney, J., and Ritchie, R. (2003). Mechanistic fracture criteria for the failure of human cortical bone. Nat. Mat. 2, 164–168. doi:10.1038/nmat832

PubMed Abstract | CrossRef Full Text | Google Scholar

Nguyen, T. T., Breeze, J., and Masouros, S. D. (2022). Penetration of energised metal fragments to porcine thoracic tissues. J. Biomechanical Eng. 144 (7), 071002. doi:10.1115/1.4053212

CrossRef Full Text | Google Scholar

Niebur, G. L., Feldstein, M. J., Yuen, J. C., Chen, T. J., and Keaveny, T. M. (2000). High-resolution finite element models with tissue strength asymmetry accurately predict failure of trabecular bone. J. biomechanics 33 (12), 1575–1583. doi:10.1016/S0021-9290(00)00149-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Öhman, C., Dall’Ara, E., Baleani, M., Jan, S. V. S., and Viceconti, M. (2008). The effects of embalming using a 4% formalin solution on the compressive mechanical properties of human cortical bone. Clin. Biomech. 23 (10), 1294–1298. doi:10.1016/j.clinbiomech.2008.07.007

CrossRef Full Text | Google Scholar

Reilly, D. T., Burstein, A. H., and Frankel, V. H. (1974). The elastic modulus for bone. J. biomechanics 7 (3), 271–275. doi:10.1016/0021-9290(74)90018-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Reilly, D. T., and Burstein, A. H. (1975). The elastic and ultimate properties of compact bone tissue. J. biomechanics 8 (6), 393–405. doi:10.1016/0021-9290(75)90075-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Saito, M., and Marumo, K. (2015). Effects of collagen crosslinking on bone material properties in health and disease. Calcif. Tissue Int. 97 (3), 242–261. doi:10.1007/s00223-015-9985-5

PubMed Abstract | CrossRef Full Text | Google Scholar

San Antonio, T., Ciaccia, M., Müller-Karger, C., and Casanova, E. (2012). Orientation of orthotropic material properties in a femur FE model: A method based on the principal stresses directions. Med. Eng. Phys. 34 (7), 914–919. doi:10.1016/j.medengphy.2011.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanborn, B., Gunnarsson, C. A., Foster, M., and Weerasooriya, T. (2016). Quantitative visualization of human cortical bone mechanical response: Studies on the anisotropic compressive response and fracture behavior as a function of loading rate. Exp. Mech. 56 (1), 81–95. doi:10.1007/s11340-015-0060-y

CrossRef Full Text | Google Scholar

Schileo, E., Taddei, F., Cristofolini, L., and Viceconti, M. (2008). Subject-specific finite element models implementing a maximum principal strain criterion are able to estimate failure risk and fracture location on human femurs tested in vitro. J. biomechanics 41 (2), 356–367. doi:10.1016/j.jbiomech.2007.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmitt, K. U., Niederer, P. F., Cronin, D. S., Morrison, B., Muser, M. H., and Walz, F. (2019). Trauma biomechanics: An introduction to injury biomechanics. 5th Ed. Cham: Springer. doi:10.1007/978-3-030-11659-0

CrossRef Full Text | Google Scholar

Sharir, A., Barak, M. M., and Shahar, R. (2008). Whole bone mechanics and mechanical testing. Veterinary J. 177 (1), 8–17. doi:10.1016/j.tvjl.2007.09.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Shore, S. W., Unnikrishnan, G. U., Hussein, A. I., and Morgan, E. F. (2013). “Bone biomechanics,” in Orthopedic biomechanics. Editor B. A. Winkelstein. 1st Ed. (Boca Raton: CRC Press), 3–48.

Google Scholar

Stefan, U., Michael, B., and Werner, S. (2010). Effects of three different preservation methods on the mechanical properties of human and bovine cortical bone. Bone 47 (6), 1048–1053. doi:10.1016/j.bone.2010.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Strømsøe, K., Høiseth, A., Alho, A., and Kok, W. L. (1995). Bending strength of the femur in relation to non-invasive bone mineral assessment. J. Biomechanics 28 (7), 857–861. doi:10.1016/0021-9290(95)95274-9

CrossRef Full Text | Google Scholar

Tang, T., Ebacher, V., Cripton, P., Guy, P., McKay, H., and Wang, R. (2015). Shear deformation and fracture of human cortical bone. Bone 71, 25–35. doi:10.1016/j.bone.2014.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Turner, C. H. (2006). Bone strength: Current concepts. Ann. N. Y. Acad. Sci. 1068 (1), 429–446. doi:10.1196/annals.1346.039

PubMed Abstract | CrossRef Full Text | Google Scholar

Untaroiu, C. D., Yue, N., and Shin, J. (2013). A finite element model of the lower limb for simulating automotive impacts. Ann. Biomed. Eng. 41 (3), 513–526. doi:10.1007/s10439-012-0687-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Ural, A., and Vashishth, D. (2006). Cohesive finite element modeling of age-related toughness loss in human cortical bone. J. biomechanics 39 (16), 2974–2982. doi:10.1016/j.jbiomech.2005.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Ural, A., Zioupos, P., Buchanan, D., and Vashishth, D. (2011). The effect of strain rate on fracture toughness of human cortical bone: A finite element study. J. Mech. Behav. Biomed. Mater. 4 (7), 1021–1032. doi:10.1016/j.jmbbm.2011.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Wolfram, U., and Schwiedrzik, J. (2016). Post-yield and failure properties of cortical bone. Bonekey Rep. 5, 829. doi:10.1038/bonekey.2016.60

PubMed Abstract | CrossRef Full Text | Google Scholar

K. H. Yang (Editor) (2017). Basic finite element method as applied to injury biomechanics (London: Elsevier). doi:10.1016/C2015-0-06702-8

CrossRef Full Text | Google Scholar

Zhai, X., Gao, J., Nie, Y., Guo, Z., Kedir, N., Claus, B., et al. (2019a). Real-time visualization of dynamic fractures in porcine bones and the loading-rate effect on their fracture toughness. J. Mech. Phys. Solids 131, 358–371. doi:10.1016/j.jmps.2019.07.010

CrossRef Full Text | Google Scholar

Zhai, X., Guo, Z., Gao, J., Kedir, N., Nie, Y., Claus, B., et al. (2019b). High-speed X-ray visualization of dynamic crack initiation and propagation in bone. Acta biomater. 90, 278–286. doi:10.1016/j.actbio.2019.03.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Zysset, P. K., Dall'Ara, E., Varga, P., and Pahr, D. H. (2013). Finite element analysis for prediction of bone strength. Bonekey Rep. 2, 386. doi:10.1038/bonekey.2013.120

PubMed Abstract | CrossRef Full Text | Google Scholar

Nomenclature

CDM Continuum Damage Mechanics

CII Crash Induced Injuries

FE Finite Element

HBMs Human Body Models

MLT Matzenmiller, Lubliner and Taylor

X-FEM Extended Finite Element Method

C₁, C₂, C₃, K_m, n_m Modified Mohr Coulomb Model fracture parameters

E_i Young’s modulus in i direction

G_ij Shear modulus in ij plane

K Bulk Modulus

S_i Ultimate Strength in i direction

S_ij Ultimate Strength in ij plane

ν_ij Poisson’s ratio in ij plane

m_ij MLT constitutive model damage rate exponent

ω_ij Damage value associated with a specific material direction

g_ij Damage rate corresponding to wij

η Triaxiality

ε_ij Current strain

ε̇_ij Current strain rate

ε_ij,f Material failure strain

ε̅_f Effective strain at failure