Computational homogenisation based extraction of transverse tensile cohesive responses of cortical bone tissue

Abstract

The numerical assessment of fracture properties of cortical bone is important in providing suggestions on patient-specific clinical treatments. We present a generic finite element modelling framework incorporating computational fracture approaches and computational homogenisation techniques. Finite element computations for statistical volume elements (SVEs) at the microscale are performed for different sizes with random osteon packing with a fixed volume fraction. These SVEs are loaded in the transverse direction under tension. The minimal SVE size in terms of ensuring a representative effective cohesive law is suggested to be 0.6 mm. Since cement lines as weak interfaces play a key role in bone fracture, the effects of their fracture properties on the effective fracture strength and toughness are investigated. The extracted effective fracture properties can be used as homogenised inputs to a discrete crack simulation at macroscopic or structural scale. The extrinsic toughening mechanisms observed in the SVE models are discussed with a comparison against experimental observations from the literature, giving beneficial insights to cortical bone failure.

Appendices

Appendix A: Hill-Mandel condition for WPBCs

This section aims to prove that the WPBCs respect the Hill-Mandel condition. According to the first-order CH, the microscale displacement field is given as

$$\begin{aligned} {\mathbf {u}} = {\mathbf {u}}_{\textsc {M}}+\tilde{{\mathbf {u}}}= \varvec{\varepsilon }_{\textsc {M}}\cdot {\mathbf {x}}+\tilde{{\mathbf {u}}}\;, \end{aligned}$$

(8)

where ${\mathbf {u}}_{\textsc {M}}$ and $\tilde{{\mathbf {u}}}$ represent the uniform deformation corresponding to the prescribed macroscale strain $\varvec{\varepsilon }_{\textsc {M}}$, and the displacement fluctuations induced by the microstructure, respectively.

The WPBCs can be regarded in an integral sense as a weak implementation of the strong PBCs that are applied pointwise. In the weak setting, the WPBCs can be expressed as (Larsson et al. 2011)

$$\begin{aligned} \int _{\Gamma ^+} \delta {\mathbf {t}}_\lambda \cdot \llbracket \tilde{{\mathbf {u}}}\rrbracket _{\Gamma } \,{\mathrm {d}}\Gamma = {\mathbf {0}} \;, \end{aligned}$$

(9)

where $\delta {\mathbf {t}}_\lambda$ can be regarded as a test function and interpreted as the weighted boundary traction vector which is unknown at the moment. The proof is given as follows:

$$\begin{aligned} \begin{aligned} \langle \varvec{\sigma } : \delta \varvec{\varepsilon } \rangle&=\frac{1}{|\varOmega |} \left( \int _{\varOmega } \varvec{\sigma } : \delta \nabla {\mathbf {u}}_{\textsc {M}}\,{\mathrm {d}}\varOmega + \int _{\varOmega } \varvec{\sigma } : \delta \nabla \tilde{{\mathbf {u}}} \,{\mathrm {d}}\varOmega \right) \\&=\frac{1}{|\varOmega |} \int _{\varOmega } \varvec{\sigma } : \delta \nabla {\mathbf {u}}_{\textsc {M}}\,{\mathrm {d}}\varOmega + \frac{1}{|\varOmega |} \int _{\Gamma } {\mathbf {t}} \cdot \delta \tilde{{\mathbf {u}}} \,{\mathrm {d}}\Gamma \\&=\varvec{\sigma }_{\textsc {M}}: \delta \varvec{\varepsilon }_{\textsc {M}}+ \frac{1}{|\varOmega |} \int _{\Gamma ^+} {\mathbf {t}} \cdot \delta \llbracket \tilde{{\mathbf {u}}} \rrbracket _{\Gamma } \,{\mathrm {d}}\Gamma \\&=\varvec{\sigma }_{\textsc {M}}: \delta \varvec{\varepsilon }_{\textsc {M}}\;, \end{aligned} \end{aligned}$$

(10)

where the equilibrium equation of SVE, the divergence theorem, the anti-periodicity of boundary tractions and Eq. (9) have been used in sequence. The WPBCs can be implemented in the context of FE by means of introducing Lagrange multipliers into the system. In the case of the multiscale crack model using an extended CH and WPBCs, the energetic equivalence across two scales has been demonstrated in (Xing and Miller 2021).

Appendix B: Comparison with a clay/epoxy nanocomposite fracture study

In the work by Msekh et al. (2018), the authors predicted fracture properties of clay/epoxy nanocomposites with interphase zones using a phase field model for fracture. The simulated nanocomposite system shares similarities to the system in this contribution. For example, inclusions (clay platelets and osteons) are stiffer than the matrix (epoxy and interstitial matrix), and both the matrix cracking and interface debonding are considered without failure in the inclusions. On the other hand, the authors used the concurrent multiscale method as compared to the semi-concurrent multiscale method herein. Some interesting findings in (Msekh et al. 2018) were (i) the overall tensile strength increased with decreasing interface thickness in the case of weaker interface compared to the matrix; (ii) The global energy release rate was significantly increased by increasing the critical fracture energy of the interface; (iii) The fracture energy of the interface did not affect the nanocomposite’s tensile strength. Meanwhile, the dissipation energy due to fracture was increased drastically with an interface rather than the fully bonded nanocomposite. Our results match their points (ii) and (iii), both of which emphasise the important role of the weak interface. However, we cannot replicate their point (i), because the interface around osteons is simulated with interface elements assuming a zero thickness.