Constitutive relationships for osteonal microcracking in human cortical bone using statistical mechanics

Abstract

The post-elastic mechanical behavior of cortical bone, which is represented by extensive microcracking once the elastic regime is exceeded, has been characterized by a nonlinear constitutive relationship for osteonal microcracking. The relationship/model is based on the formalism of Statistical Mechanics, allowing the degree of irreversibility to be calculated using the increase in entropy associated with the progression of microcracking. Specific tensile and bending tests were conducted to compare theoretical predictions of constitutive relationships to empirical curves. In addition, the tests were utilized to determine the model’s parameters, whose values were used to explicitly calculate the entropy increase. A large sample was used: 51 cortical bone coupons (dog-bone-shaped specimens) were extracted from the 4th ribs of numerous individuals and subjected to uniaxial tensile testing. Additionally, fifteen complete 4th ribs were used for bending tests. Displacement and strain fields were measured for both types of tests using digital image correlation or video recordings of the tests. All experimental specimen data were successfully fitted to the model, and all constitutive parameter values were found to be correlated with anthropometric variables. Explicit entropy calculations indicate that microcracking is minimal for low strain and, initially, stress is nearly proportional to strain. After a certain point, significant microcracking occurs, and the relationship between stress and strain becomes invalid. Several significant associations between constitutive parameters and age have also been identified.

1 Introduction

The process of microcracking in cortical bone is a thermodynamically irreversible process that eventually leads to macroscopic fracture of the bone. Such cracking depends on the stress and strain inside the bone. In biological tissues, stress depends on mechanical loads which, at the material level, are mathematically related to strain by means of a constitutive material model. Thus, when a load is imposed on a bone, the material reacts with a stress–strain response. This response is intrinsic to the material and depends on the magnitude, the rate, the direction and the duration of the applied load.

A number of comprehensive constitutive mechanical models have been published to mathematically relate stress and strain. Some of these models use plasticity theory [1,2,3,4,5,6,7], continuous damage theory [8,9,10,11,12,13], rheological theories [14, 15] (including nonlinear viscoelastic and viscoplastic models [16,17,18]), microstructural parameters and microcontinuum models [19,20,21], statistical models such as the random spring network model (RSNM) [22, 23] and even, a combination of the above characteristics [16, 24,25,26,27].

Although the vast majority of mechanical constitutive models have obvious practical value, only a few of them analyze the thermodynamic aspects of the cracking process [28, 29]. Specifically, few of the proposed constitutive models allow an explicit calculation of the entropy associated with cracking, which, as an irreversible process, must increase as microcracking progresses within the bone. Despite the fact that a certain number of models describe the microstructure of cortical bone [30, 31], the progression of microcracking has received less attention than other aspects influencing mechanical behavior [32]. In particular, the microdamage accumulation and the presence of microcracks have received some attention [33, 34], as well as the occurrence of plastic strain [31, 35, 36]. However, not all published constitutive models address the microcracking progression and sometimes, models are not accompanied by a great amount of experimental evidence to validate its performance [16].

The purpose of this study is to develop a thermodynamically consistent mechanical model of bone microcracking progression. The model’s formulation allows for the explicit calculation of the entropy increase associated with cracking progression (verifying that it is an irreversible process). The innovative aspect of this paper is the use of Statistical Mechanics to model the microcracking process and calculate the entropy associated with the irreversibility of the damage. Moreover, in this study the model is validated with experimental data from tensile and bending tests of human rib cortical bones, computing the parameters of the model and determining its relation with anthropometric variables, as well as calculating the entropy associated with the microcracking that evolves in the microstructure during the tests.

2 Data methods and proposed model

This section describes the theoretical basis for the constitutive relationships as well as the experiments (tensile and bending tests) used to validate their applicability. The model was used to compute the stress–strain relationship and the entropy increase for a monotonic mechanical load application. The results were compared to experimental data. The model does not account for bone remodeling because it assumes that stress will increase until bone fracture occurs.

2.1 Microcracking and statistical mechanics

For the purpose of deriving the microcracking constitutive relationship, it is assumed that bone is composed of osteons, which serve as the model’s fundamental elements (for this reason, lower hierarchical levels will not be analyzed). Osteons are aligned predominantly along the longitudinal direction of long bones, resulting in differential mechanical behavior depending on direction [37,38,39,40]. Thus, cortical bone can be considered as a transversely isotropic material.

To capture the behavior of the nonelastic region, a transversely isotropic mechanical model was developed for Strain Energy Density Functions (SEDF) using the Canonical Ensemble Formalism of Statistical Mechanics. In general, the space between osteons is weaker than the osteons themselves, so there are more microcracks in the inter-osteon space [41,42,43]. When cortical bone is subjected to mechanical forces in the directions of dominant tensile stress, this must be accounted for in a constitutive model and even a fracture model [44, 45].

The Canonical Ensemble Formalism permits for the introduction of the quenched-disorder parameter through the partition function $Z(\beta ,\textbf{E})$ [46,47,48]. This quenched-disorder parameter describes the initial microcracking in the specimen (prior to stress application), and the constitutive relationships allow measuring the progression of microcracking by computing the entropy. Microcracking is modeled as a propagation in an osteon lattice for this purpose. Therefore, the model relates the energy introduced into the material to the propagation of microcracks within the osteon-formed microstructure. This microcracking is analyzed using combinatorics, which computes the number of potential paths along the lattice of osteons.

2.1.1 Canonical ensemble and Helmholtz free energy

The collection of osteons constitutes an ensemble that forms the cortical bone tissue, and this ensemble is amenable to a statistical mechanical description. In the Statistical Mechanics formalism, the starting point is a large number of elementary components that interact by exchanging energy and influencing each other. Several macroscopic magnitudes characterize a system’s macrostate (in our case, will be the strain and the total microcracking length), and each macrostate is compatible with many microstates. Based on the number of microstates compatible with each macrostate, Statistical Mechanics attempts to determine the most probable macrostate, which is proportional to the number of compatible microstates, and the average behavior for each macrostate [49,50,51].

Before loading the specimen, the initial state is a partially cracked macrostate with low microcracking, irregularities and some initial defects. As the load increases, so does the strain, and consequently, the mechanical energy and microcracking in the bone increases. The cracks propagate primarily in the inter-osteon space, through the interstitial bone and cement lines [52]. To calculate the total microcracking length (TML), a plane orthogonal to the surfaces which contain cracks is considered (this plane also contains the barycentric axis of the rib bone). The intersection of all microcracks with such a plane yields the TML (adding lengths for all the cracks contained in the plane). This TML ($\ell$) can be compared to the average length of the osteons $\ell _0$, and the cracking number $k = \ell /\ell _0$ is then calculated. The cracking number identifies the cracking macrostate for each strain level ($\varepsilon$), having an energy of $E_k = E(k,\varepsilon )$. For a given cracking number k, different configurations or paths are possible for each crack in the osteon lattice.

To determine the partition function $Z(\beta ,\textbf{E})$, it is necessary to count the number of potential cracking configurations in the osteon lattice [46, 50, 51]. The term cracking microstate will refer to any microscopic configuration that is compatible with a given macrostate or strain. As a result, counting the number of microstates entails counting the number of configurations that are possible depending on the cracking number k (or, alternatively, the TML $\ell = k\ell _0$). The partition function $Z(\beta ,\textbf{E})$ is used to calculate the probability $p_k$ of each macrostate with cracking number k, using the traditional Maxwell–Boltzmann statistics; specifically [50, 51]:

$$\begin{aligned} p_k=\frac{g_k\ e^{-\beta E_k}}{Z(\beta ,\textbf{E})},\qquad Z(\beta ,\textbf{E})=\sum _{j=1}^{\Omega }e^{-\beta E_j(\textbf{E})} \end{aligned}$$

(1)

where $\Omega$ is the total number of microstates. Then, the total energy of the system $E = \sum _j p_jE_j$ is given by the strain field represented by the tensor $\textbf{E}$, and $\beta$ is the quenched-disorder parameter of the system [47, 48]. In the lattice constituted by osteons and their junctions, the increase in microcracking will lead to an increase in disorder and, therefore, an increase in entropy. However, unlike thermomechanical systems where $\beta =1/T$ the increase in entropy here is not associated with any change in temperature, but with the progression of other dissipative and irreversible processes. The partition function allows obtaining the strain energy function that coincides with the free Helmholtz energy density:

$$\begin{aligned} \Psi (\beta ,\textbf{E}) = -\frac{1}{\beta }\ln Z(\beta ,\textbf{E}) \end{aligned}$$

(2)

This function coincides with the SEDF, and from it, all the necessary constitutive relations can be derived.

2.1.2 Microstates, macrostates and energy levels

As previously stated, each load level carries a well-defined strain which is intrinsic to the corresponding macrostate’s specific energy $E_k$. In addition, for each macrostate, the number of possible crack configurations (microstates) in the osteon lattice will be determined. Consequently, $Z(\beta ,\textbf{E})$ is defined by a sum over all possible microconfigurations. As in most systems, the number of microstates increases as the total energy and, consequently, the strain level increases. If the number of microstates or configurations compatible with each energy level is known, however, Z can be rewritten as a sum over the various energy levels:

$$\begin{aligned} Z(\beta ,\textbf{E})=\sum _{k} g(E_k)e^{-\beta E_k (\textbf{E})} \end{aligned}$$

(3)

where $g(E_k)$ is a function that counts the number of cracking paths compatible with the energy level $E_k$ (which will be studied in the next section). Then, it follows that as the load increases, both the energy $E_k$ and the TML increase, since each individual crack in the osteon lattice can take different paths. Thus, each macrostate k specified by $E_k$ implies a new stage of crack propagation that lengthens it. We assume that the energy of each macrostate in expression (3) is given by $E_k=k\psi _0$, where k is the cracking number and $\psi _0$ is the SEDF of a transversely isotropic material (without cracking), calculable as a combination of the following invariants [53, 54]:

$$\begin{aligned} \psi _0(\textbf{E}) = \psi _0(I_1(\textbf{E}), I_2(\textbf{E}), I_3(\textbf{E}); I_4(\textbf{E},\varvec{a}), I_5(\textbf{E},\varvec{a})) \end{aligned}$$

(4)

where $\varvec{a}$ is a vector field giving the preferential alignment direction of the osteons along the rib, and the invariants $I_i$ are the trace or linear invariant ($I_1$), the quadratic invariant ($I_2$), the determinant ($I_3$), and the linear and quadratic invariants representing transversal anisotropy ($I_4,I_5$):

$$\begin{aligned} \begin{array}{rl} I_1(\textbf{E}) &{}= \text {tr}(\textbf{E})\\[1ex] I_2(\textbf{E}) &{}= \frac{1}{2}[\text {tr}^2(\textbf{E})-\text {tr}(\textbf{E}^2)],\\[1ex] I_3(\textbf{E}) &{}= \det (\textbf{E}) \end{array}\quad \begin{array}{rl} I_4(\textbf{E},\varvec{a}) &{}= \varvec{a}\cdot \textbf{E}\varvec{a}\\[1ex] I_5(\textbf{E},\varvec{a}) &{}= \varvec{a}\cdot \textbf{E}^2\varvec{a} \end{array} \end{aligned}$$

(5)

Given the low strain values for cortical bone, any function of strain can be adequately approximated by a polynomial of small degree. Because the highest order invariant $I_3$ is cubic in strain tensor in this case, a polynomial of degree three which is a combination of all the above invariants will be considered (as it is the smallest degree that permits the use of all invariants):

$$\begin{aligned} \psi _0\approx \mu _2I_1^2+\mu _3I_1^3+\nu _2I_2+\nu _3I_1I_2+\varsigma _3I_3+\rho _2I_4^2+\rho _3I_4^3+\tau _2I_5+\tau _3I_4I_5 \end{aligned}$$

(6)

where terms with degree higher than three in strain have not been considered due to the limited low values of strain in bone, and $\mu _i, \nu _i, \rho _i, \tau _i$ are parameters of the model that must be fitted. The parameters with subscript 2 ($\square _2$) multiply quadratic in the strain components; those of subscript 3 ($\square _3$) multiply cubic in the strain. Note that, unlike $\Psi$ in Eq. (2), $\psi _0$ does not incorporate the progressive degradation of the specimen.

It is also worth noting that neither internal damage variables $\varvec{\alpha }$ nor evolution equations $\varvec{\dot{\alpha }} = f(\varvec{\alpha },\textbf{E})$ are specifically included [55]. This is due to the fact that, in a mechanical test with monotone loading at constant velocity, as in our case, the equation of damage evolution is integrable and can be expressed only in terms of longitudinal deformation $\varepsilon$. Use the chain rule for time derivative to see it:

$$\begin{aligned} \frac{\text {d}\varvec{\alpha }}{\text {d}t} = \frac{\text {d}\varvec{\alpha }}{\text {d}\varepsilon }\dot{\varepsilon } = f(\varvec{\alpha },\textbf{E}(\varepsilon )) \end{aligned}$$

(7)

and $\textbf{E}(\varepsilon )$ because of relationships (23) and (24). Integrating this last equation we have $\varvec{\alpha } = \Phi (\varepsilon ) = \hat{\Phi }(\textbf{E})$. In other words, we are not proposing a comprehensive constitutive model here [32, 56], but rather a constitutive relationship that can be compared to the tensile and bending tests we have conducted.

2.1.3 Counting microstates: possible paths for microcracks

With the above considerations, partition function (3) is constructed, which in turn allows computing the SEDF $\Psi$ using expression (2), which finally gives the constitutive relation sought. This procedure requires first calculating the numbers $g_k = g(E_k)$ that give the number of microstates compatible with the energy $E_k$ of each macrostate. The number $g_k$ is given by the number of possible paths for a crack in osteon lattice (number of microstates compatibles with $E_k$ and $\ell = k\ell _0$).

Considering a lattice with osteons, separated by bonds corresponding to the osteon interface, microcracks are considered to propagate along these interfaces, in a way that a TML $\ell = k\ell _0$ developed in k steps.

Figure 1 compares three regular lattices in order to calculate the number of paths $g_k$ for a specific cracking number k or energy $E_k$. The first lattice is an orthogonal lattice, whose cells represent osteons and lines represent interfaces between them. For each phase of crack advancement, there are three possible paths (Fig. 1a). The total number of paths is $3\cdot 3\cdot 3\dots =3^m$, where m is the total number of steps of propagation. However, cracks in a material have a tendency to advance and not to make abrupt direction changes that result in a higher energy level. Consequently, limiting the retreat of the crack advance (Fig. 1b), there are $2\cdot 2\cdot 2\dots =2^m$ possible paths (microconfigurations) compatible with a given TML $\ell$.

Analyzing a triangular lattice, the possible paths are $5\cdot 5\cdot 5\dots =5^m$, whereas these paths are $3\cdot 3\cdot 3\dots =3^m$ without backward movements (Fig. 1c, d). In a hexagonal lattice, there are $2\cdot 2\cdot 2\dots =2^m$ paths in the general case, but $2\cdot 1\cdot 2\cdot 1\cdot 2\dots =(\sqrt{2})^m$ in the case for a crack that only advances (Fig. 1e, f).

Fig. 1

Three idealized lattices: Orthogonal: a the crack encounters three possible paths $g_k=3^k$ b without backward movements; the crack has two permissible directions such that $g_k=2^k$. Triangular: c the crack has five possible paths $g_k=5^k$, d the crack has three permissible directions without backward movements, and $g_k=3^k$. Hexagonal: e the crack has two possible paths $g_k=2^k$, f for no backward movement; there are $2\cdot 1\cdot 2\dots$ possible paths and $g_k=(\sqrt{2})^{k}$

According to the previous analysis, the number of possible paths is of the form $g_k=r^k$, where r is the lattice parameter. In the lattices analyzed, r satisfies $\sqrt{2}\le r \le 3$. The maximum number of advance steps is defined as N (assuming the average cortical thickness is 0.8 mm and the diameter of osteons is 0.2 mm, $N \approx 5$). However, a lattice of osteons is less ideal than regularly analyzed lattices, so it will be established that 1$< r \le$ 3 and r and N will be determined from the fitting procedure, see Sect. 2.3.2.

2.1.4 Proposed model for microcracking

After obtaining the function $g_k=r^k$ that defines the number of possible microcrack paths, partition function (3) is given by:

$$\begin{aligned} Z(\beta ,\textbf{E}) = \sum _{k=1}^{N} r^ke^{-\beta k\psi _0(\textbf{E})} = re^{-\beta \psi _0}\ \frac{1-(re^{-\beta \psi _0})^N}{1-re^{-\beta \psi _0}} \end{aligned}$$

(8)

where the sum goes up to N, being N the total number of steps of crack advance (which cannot be infinite), and $E_k=k\psi _0$. Clearly, Z can then be rewritten as a closed-form expression. Then, the SEDF $\Psi$ is obtained by using (4) and deriving $\Psi = (1/\beta )\ln Z$. Thus, the stress tensor can be calculated as follows:

$$\begin{aligned} \textbf{S}(\textbf{E}) =\frac{\partial \Psi (\beta ,\textbf{E})}{\partial \textbf{E}}= -\frac{1}{\psi _0}\frac{\partial \ln Z}{\partial \beta }\ \frac{\partial \psi _0}{\partial \textbf{E}} \end{aligned}$$

(9)

from Eqs. (2) and (8). Additional details related to the convenient change of derivation variable $\textbf{E}$ to $\beta$ are described in appendix section 5.1. Using the same details of the appendix, it can be seen that the constitutive relationships take the form:

$$\begin{aligned} \textbf{S}(\textbf{E}) = \bar{\textbf{S}}(\textbf{E})\left( \frac{1}{1-re^{-\beta \psi _0}}- \frac{Nr^Ne^{-\beta \psi _0N}}{1-r^Ne^{-\beta \psi _0N}}\right) = \bar{\textbf{S}}(\textbf{E})\cdot \Phi (\beta ,\psi _0(\textbf{E}),N,r) \end{aligned}$$

(10)

where $\bar{\textbf{S}}$ is the “crackless stress,” see Eq. (22), and $\Phi$ is the crack function that corrects $\bar{\textbf{S}}$ to account for microcracking effects. Moreover, $\psi _0$ is the “crackless SEDF” defined in (6), $\beta$ is the quenched-disorder parameter, and N, r are parameters that define the possible paths that a crack can take in the osteon lattice; specifically, r is the lattice parameter ($1<r \le 3$) and N is of the same order as the number of osteon layers ($N\approx 5$). Therefore, the parameters to be determined are $\varvec{\theta } = (\mu _2,\mu _3,\nu _2,\nu _3,\varsigma _3,\rho _2,\rho _3,\tau _2,\tau _3)$, which are included in the strain energy function $\psi _0$ along with the additional parameters $\beta , N, r$. Once these parameters are fitted, the Canonical Ensemble Formalism allows computing the derivative of SEDF (2):

$$\begin{aligned} H(\textbf{E},\beta ) = \beta ^2 \frac{\partial \Psi }{\partial \beta } = \ln Z(\textbf{E}_f,\beta ) - \beta \frac{\partial \ln Z(\textbf{E}_f,\beta )}{\partial \beta } \end{aligned}$$

(11)

In Sect. 3 explicit computations of the increment in entropy will be presented.

2.2 Materials and experimental setting

Author-obtained experimental data were used to evaluate the suitability of the proposed model. The experimental data are derived from both tensile tests of small specimens of cortical bone and bending tests of complete ribs. These tests allowed the model’s parameters to be fitted, which enables the calculation of several additional quantities, including entropy. The subsequent sections describe the material used, its origin, the specifics of the mechanical tests and the mathematical procedure of numerical fitting.

Specimens of human ribs were obtained from forensic autopsies conducted at the Forensic Pathology Service of the Legal Medicine and Forensic Science Institute of Catalonia (IMLCFC). According to the protocol approved by the monitoring committee of the collaboration project between our research center and the forensic institute, all specimens were obtained from individuals whose thorax-related deaths were not the result of trauma.

Tensile and bending tests were performed with a MicroTEST EM2/20 (MicroTEST^®, Madrid), and the force was measured with a 500-N HBM S9M load cell attached to a Spider 8–30 acquisition system of HBM^® (Spider, Darmstadt). The strain rate was in both cases less than $\dot{\varepsilon } = 0.02$ sec$^{-1}$ to reproduce nearly quasistatic conditions. The tests were recorded with a high-speed (60 Hz) and high-resolution (megapixel) camera (PCO 1200s) for the subsequent processing of the images and the calculation of strain.

The displacement and strains were calculated using changes in the position of the specimen material points from one image of the video to the next, allowing the creation of extremely accurate curves with minimal experimental noise and the calculation of displacement and strain fields with high precision, using digital image correlation (DIC), as described in Sect. 2.3.1.

2.3 Mechanical testing

In this section, the details of the mechanical tests, the number of specimens, and the experimental set used in each case are provided:

• Uniaxial tensile tests. For the tensile tests, $n_T=$ 51 specimens (35 males and 16 females) aged 55±20 years (from 21 to 91 years) were used. According to the procedure described in [57,58,59], these ribs were machined to obtain coupons from the anterior rib region. First, soft tissue was removed and rib cortical slices were obtained with a low-speed diamond saw from the anterior portions of the rib. Then, two holes of 2 mm were performed on the slices, at a distance of 20 mm between each other (for centering of the slice and fixation in the clamping system of the tensile test). Lastly, coupon shape was machined with a milling machine and meticulously polished to ensure a constant gauge thickness. The specimens measured 28 mm in length and 8 mm in width, with a gauge length of 6 mm and a gauge width of 2 mm and an approximate thickness of 0.5 mm (Fig. 2a). This size is small enough to ensure a constant cross section in all the gauge length of the specimen. In the tensile tests, the coupons were placed in the machine and the deformation was increased until the specimen was broken by the reference length. From the set of images of the video recorded during the test, the displacements of various points on the specimen throughout the test were obtained through a DIC procedure of Matlab^®, according to the methodology of [59]. To do this, a rectangular mesh of points was created in the reference length and the displacements of these points were determined during the test (previously, a random dot pattern was painted on the specimen to identify the position of these points), see Fig. 2b.

• Plane Bending tests. Each rib was subjected to bending within a plane, using a methodology developed by the authors previously [60]. Before the tests, soft tissue and costal cartilage were removed. The complete ribs were placed in the plane of the machine, introducing their ends in a U-guide placed on the upper bench and the outer central region in contact with the actuator that exerts the force. The upper guide was covered with lubricant, and the rib ends were wrapped in polytetrafluoroethylene tape to minimize friction. Thus, the ends of the rib could slide freely in the guide while it was pressed on the center of the rib, causing a bending moment that produced a loss of curvature of the rib, which resulted in the appearance of a traction zone and a compression zone. Moreover, each rib was placed between four safety rods, which have not been in contact with the specimen. For the tests $n_B=$ 15 complete ribs (11 men and 4 women) from a sample ranging in age from 26 to 62 years have been utilized (sample average 51±11 y. o).

As described in Sect. 2.3.1, both tests were captured on video and the frames were analyzed using DIC methods, which enabled the placement of a series of reference points throughout the trial in order to determine the displacement field and strain field.

2.3.1 Data processing and digital image correlation

This section describes the strain calculation for tensile and bending tests. It is possible in both instances to determine stress and strain from two-dimensional images due to the absence of stress in the direction perpendicular to the recording plane (explicit forms of the strain tensor are given in section 5.2). In both situations, the position of a mesh of material points is explicitly measured. The displacement field is computed using interpolation from a collection of measured point coordinates. By comparing the obtained displacements to the initial coordinates of the points, it is possible to compute the deformation gradient and strain field numerically.

Fig. 2

a Coupon dimensions for the specimens machined from rib cortical bone. b The initial grid used for defining a collection of material points which will be tracked during the elongation of the specimen. Points are separated between them 50 pixels

As illustrated in Fig. 2, the uniaxial tensile tests use a mesh of orthogonal points. A DIC motion tracking procedure identifies specific patterns around each point, and the location of each material point of the mesh is determined for each instant of time (60 readings are made per second and 200 points are used, so it is easy to interpolate in time and space with a high degree of precision). From the instantaneous point positions, the displacements $\varvec{u}(\varvec{x})$ were determined, and the finite deformation theory was applied to calculate the components of the Green–Lagrangian strain tensor as follows:

$$\begin{aligned} E_{ij} = \frac{1}{2}\left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_i}{\partial x_j} + \sum _k \frac{\partial u_k}{\partial x_i}\frac{\partial u_k}{\partial x_j} \right) \end{aligned}$$

(12)

for $i, j \in \{1,2\}$ [being $x_1 = x$ and $x_2 = y$]. In this way, the longitudinal and transverse deformations ($E_{xx} = E_{\text {long}}$ and $E_{yy}=E_{zz} = E_{\text {trans}}$, respectively), related to the applied stress level, were determined. The stress was calculated at each time with Eq. (10) (it should be noted that even when we compute the real or instant area, it is always close to the initial area, because cortical tissue strain is very low, typically 3%).

Fig. 3

Experimental setting for the bending tests: the rib ends are inserted into the lubricated guide and the force is exerted in the lower center area. With increasing force, the ends open freely. The four safety rods protect against slipping

Figure 3 shows the experimental setting and the set of points used in the DIC procedure that gives the point positions along the outline of the specimen. As before, motion tracking identifies the points in each frame as the rib is deformed (the XY plane coincides with the frame plane). Thus, the displacements of various points on the upper and lower contour of the specimen were determined and the displacements of the midline y(x) between these contours were computed, which was considered as the barycentric line of the rib (it was verified with a CT that the error committed between considering the barycentric line real and the median line was $3\pm 2.5\%$ in the central zone). From the displacements and the current midline $y_t(x)$, the current curvature $\chi _t$ of the midline was computed as:

$$\begin{aligned} \chi _t = \frac{y''_t(x)}{[1+{y'_t}^2(x)]^{3/2}} \end{aligned}$$

(13)

from this curvature, the strain $\varepsilon _t$ was calculated using Eqs. (25) and (26). For this computation the Frenet-Serret coordinates $(\xi ,\eta ,\zeta )$ were used, following the procedure described in a previous work [60]. The midline is the set of points with $\xi = 0$, while $\eta$ is an orthogonal coordinate to midline (following the direction of the normal vector $\varvec{n}$ to midline), and $\zeta$ is the coordinate orthogonal to $\xi$- and $\eta$-coordinate lines (thus, $\zeta$ grows in the direction of the binormal vector $\varvec{b}$ to midline), see Fig. 8.

In section 5.2 of appendix, further details are given for the computation of coordinates and the strain tensor for tensile and bending tests.

2.3.2 Fitting procedure

The force and strain were computed for each instant of time. These data are used to fit the model’s constituent parameters: $\varvec{\theta } = (\mu _i,\nu _i,\rho _i,\tau _i,\varsigma _i)$ which appear in Eq. (6), along with the $(r,N,\beta )$ which appear in partition function (8). For the least squares problem, some penalty functions $\Theta _K(\varvec{\theta };\beta ,r,N)$ are defined as a sum of squares, which must be minimized numerically. In the case of tensile tests, the problem is somewhat simpler, since it is possible to calculate the stress directly:

$$\begin{aligned} \Theta _T(\varvec{\theta };\beta ,r,N) = \sum _{t=1}^m \left[ \bar{S}_{xx}^{(t)} - S_{xx}(\bar{E}_{xx}^{(t)},\bar{E}_{yy}^{(t)};\varvec{\theta },\beta ,r,N) \right] ^2 \end{aligned}$$

(14)

for minimization, about $m \approx 3500$ different instants of time (this number depends on the specimen) were used. For each instant, the experimental stress $\bar{S}_{xx}^{(t)} = \sigma _t$ is computed. In addition, using a DIC procedure (see Sect. 2.3.1), the strains on the plane $\bar{E}_{xx}^{(t)}$ and $\bar{E}_{yy}^{(t)}$ were computed for each instant t; then constitutive relationship (10) allows obtaining $S_{xx}(\bar{E}_{xx}^{(t)},\bar{E}_{yy}^{(t)};\varvec{\theta },\beta ,r,N)$.

The fitting is somewhat different for tensile tests than for bending tests, because in the former it is possible to calculate the tension directly from the force and area of the cross section. For bending tests, because it is not possible to use a Navier-like formula for a nonlinear material [60], a different penalty function was used, which involved the bending moment $M_b$ associated to the pushing force:

$$\begin{aligned} \Theta _B(\varvec{\theta };\beta ,r,N) = \sum _{j=1}^N \left[ \bar{M}_b^{(t)}-M_b(\chi _t;\varvec{\theta },\beta ,r,N) \right] ^2 \end{aligned}$$

(15)

where the bending moment $M_b(\chi _t;\varvec{\theta },\beta ,r,N)$ is a expression depending on the current curvature of the rib, which can be computed by means of expression (27) in appendix Sect. 5.3.

3 Results

The suitability of the model was assessed by comparing the experimental curves from the mechanical tests with the curves obtained by fitting the parameters of constitutive relationship (10). Stress and strain data were obtained for tensile tests of $n_T=51$ coupons (dog-shaped specimens) of human rib cortical bone and for $n_B = 15$ entire 4th ribs. Five typical stress–strain curves for tensile and bending tests and a comparison with the model fittings are presented in Fig. 4. It can be clearly seen that specimens exhibit a typical nonlinear behavior, consisting in a loss of stiffness by strain increase, which is captured by the model in both tensile and bending situations.

As it is shown in the same figure, the proposed stress–strain relationship was able to reproduce properly the experimental stress–strain curves with high precision. Indeed, for all fittings, a correlation $R^2 > 0.9999$ was obtained. The stress–strain data were used to compute the constitutive parameters $\mu _i, \nu _i, \gamma _i, \varsigma _i; \beta , N, r$, for each specimen. Table 1 shows the average parameter values for tensile and bending tests. Moreover, the mean values for both test types, divided in age groups, can be seen in Tables 2 and 3 of Appendix section 5.4. In almost all instances, the estimation error for the parameters in bending tests is greater than that computed error for tensile tests. This may be partially attributable to the size-dependent effects observed for the mechanical properties of cortical bone [61], or to the limitations of bending tests [60].

Fig. 4

a Tensile tests: stress–strain experimental data (symbols) for three specimens and fittings to the model (lines). b Bending tests: stress–strain curves for two specimens (dotted lines) and fittings to the model (continuous line)

Table 1 Parameter values for the coupon specimens under uniaxial tensile tests fitted using the model proposed, in age groups

In most cases, $\mu _2, \rho _2, \beta , N$ and r had positive values, whereas the remaining parameters had negative values. In addition, r and N always took a value of $N\approx 3$ and $r\approx 1.30$, respectively.

There is a clear trend of $\beta$ to increase with age (see Table 2). This trend was supported by a linear regression analysis, which revealed a statistically significant rise in $\beta$ with age ($p<0.0001$).

In order to comprehend the relationships between the constitutive parameters, a Principal Component Analysis was conducted to examine the joint variability of the model’s constants in the sample [62]. Only four Principal Components (PC) are required to account for the 80.5% of the observed variability in the sample, with PC$_1$ representing a 38.8% and PC$_2$ a 23.9%. As seen in Fig. 5b, the largest contributions to PC$_1$ are $\mu _2$ (18.2%) and $\nu _2$ (18.5%), whereas $\rho _2$ (34.8%) and $\tau _2$ (33.7%) are the largest contributors to PC$_2$.

The main influence of $\mu _2$ and $\nu _2$ shows that the coefficients of the first two invariants $I_1(\textbf{E})$ and $I_2(\textbf{E})$, defined in Eq. (5), have a greater explanatory weight in the mechanical behavior. On the other hand, the major influence of $\rho _2$ and $\tau _2$ denotes that these coefficients of the invariants $I_4(\textbf{E},\varvec{a})$ and $I_5(\textbf{E},\varvec{a})$ that, in the case of uniaxial stress, due to the particular characteristics of the application of the stress both multiply the axial deformation, are significant in the mechanical response.

In addition, to estimate the irreversibility of the microcracking process, a computation of the entropy change was added and is displayed in Fig. 6a. As soon as the stress–strain curve is slanted downward, the model’s computed entropy increases dramatically. This appears to confirm the hypothesis that softening due to loss of stiffness is associated with irreversible changes within the cortical bone, as the amount of softening correlates clearly with changes in entropy, as shown in the graph.

Figure 6b depicts the change in entropy between $\beta$ and a generic strain measurement $\epsilon$. An increase in the quenched-disorder parameter and/or the strain level results in an increment in entropy (the increase in entropy has been calculated for each specimen as $\Delta H = H(\textbf{E}_f,\beta )-H(\textbf{0},\beta )$).

Fig. 5

a PCA of the set of parameters obtained, where CP1 represents the 38.8% of the variability and CP2 the 23.9%, and $\mu _2$ and $\nu _2$ are mainly aligned with CP1, while $\rho _2$ and $\tau _2$ are mostly aligned with CP2, b comparison for the values of microcracking parameters $(\beta , N, r)$ obtained in tensile (T) and bending (B) tests

Fig. 6

a Model fitting for one specimen and entropy increase $\Delta H$ computed from the model [light gray line, convex curve]; a fast-growing curve for the entropy is seen once the stress–strain curve [black line] loses the linear trend, and microcracking and disorder begin. b Representation of the variation of entropy with $\beta$ and a generic deformation $\epsilon$

In order to interpret these results, a regression analysis was performed with the model’s constant values. For this purpose, the actual work to fracture W was compared to the ideal-elastic work to fracture (without softening) $W_e = YE_{xx}^2/2$, where Y is the Young modulus. This allows the ratio $w = W/W_e$ to be defined. In the absence of microcracking, $w \approx 1$ is anticipated. This magnitude is significantly correlated with the number of layers of osteons N (p-value $< 0.001$), the shape factor of the osteon lattice r (p-value $< 0.0001$), the quenched-disorder parameter $\beta$ (p-value $< 0.01$) and the entropy increase $\Delta H$ (p-value $< 0.0001$). In addition, there is a significant correlation between the increase in entropy $\Delta H$ and the quenched-disorder parameter $\beta$ (p-value $\approx 0.0 1$). Other relationships between the constitutive parameters have been identified and will be discussed in the next section.

4 Discussion

A novel microcracking model of cortical bone, based on Statistical Mechanics, has been presented to predict how the progression of the microcracking process affects the stress–strain behavior (softening behavior). Given the distribution of osteons, the model assumes that the cortical bone as a whole is a statistical ensemble of osteons that can be modeled as a transversely isotropic material. The model provides an excellent fit to the empirical stress–strain curves, allowing for the computation of the degree of irreversibility with the progression of microcracking via entropy increase.

Unlike many other published theoretical models, the proposed model has been applied to a specific experimental data obtained by the authors. A large sample of cortical bone specimens is used in the experimental testings: $n_T = 51$ were tested under tensile conditions and $n_B = 15$ which were tested under bending conditions. The experimental results confirmed previous findings on the nonlinear behavior of cortical bone under axial tensile situations and more complex situations [23, 38, 42, 58]. A quasi-linear elastic behavior is found for small stress, and a loss of stiffness (strain softening) is observed after a certain level of strain. In the model, this strain-induced loss of stiffness is predicted to be associated with internal microcracking process.

The fitting of the constitutive parameters of the model to experimental data from human rib cortical bone has been used to examine the suitability of the model. The constitutive parameters obtained for both types of tests were similar, see Fig. 5b.

On the other hand, since microcracking is a dissipative and irreversible process, it is expected to lead to an increase in entropy. The fitted values of the parameters were used to calculate the entropy increment, by means of Eq. (11). In the experimental curves of Fig. 4, it can be seen that when the stress–strain behavior is approximately linear, the microcracking is practically null, and only when microcracking becomes important it is observed that the stress–strain curve becomes concave ($\text {d}^2\sigma /\text {d}\varepsilon ^2 < 0$), and thus the entropy increases. This confirm that, when the stress–strain curve shows the beginning of a marked loss of stiffness and a loss of linearity, entropy simultaneously increases rapidly. This fact suggests that the loss of stiffness is caused by microcracking (as measured by the increase in entropy $\Delta H$); therefore, the model predicts a clear correlation of both phenomena. The prediction of this correlation suggests that the approach based on Statistical Mechanics formulation is promising: it allowed to calculate numerically the entropy, while predicting accurately the stress–strain behavior and showing the interrelation between the two phenomena. In the previous section, it was mentioned that a significant correlation (p-value $< 0.0001$) was found between the ratio of work to fracture with and without softening $0< w < 1$ which directly accounts for the existence of softening, and the value of entropy changes $\Delta H$. This correlation is what allows us to state that the model adequately represents the softening due to microcracking and in a manner consistent with the thermodynamical framework. In addition, an adequate good fitting between all the experimental curves and the modeled curves is observed, thus accomplishing one of the main objectives of this research: to propose a nonphenomenological model, which can explain macroscopic stress–strain curves from a reasonable account of the microstructure modeling. Unlike other previous microcracking models, the proposed model accounts for finite strain, tissue anisotropy and thermodynamic irreversibility.

Fig. 7

Scatterplot of the quenched-disorder parameter $\beta$, related to preexisting microcracking in cortical bone tissue, against age

A clear trend in the variation of the quenched-disorder parameter $\beta$ is observed, showing a significant increase of its value with age, see Fig. 7. This suggests that the initial quenched-disorder parameter for elderly people is associated with a greater contribution of more initial states than in young people, which could explain the lower yield stress and work to fracture values (strain energy density) observed in some studies [58]. This implies that the $\Phi$, which is the term in parentheses in Eq. (10), is decreasing in $\beta$, i.e.,

$$\begin{aligned} \frac{\partial \Phi }{\partial \beta } = \frac{\partial }{\partial \beta }\left( \frac{1}{1-re^{-\beta \psi _0}}-\frac{Nr^Ne^{-\beta \psi _0N}}{1-r^Ne^{-\beta \psi _0N}}\right) < 0 \end{aligned}$$

(16)

so an increase in $\beta$ entails a lower stress level, given a strain level and other constants being equal. This prediction is reflected in the experimental data that show the expected trend for elderly people. This increase in the initial microcracking or disorder parameter agrees with the direct micrographical observations of other authors that showed that an increase in microcracking is responsible for lower mechanical properties, even in weaker bones the crack path is less deflected [43, 63]. Note that the term $\Phi$ in Eq. (10) matches the expected value of the TML $= k\ell _0$. To see this, use the probability distribution $p_k$ given by (1) for the cracking number k, then we have:

$$\begin{aligned} \overline{\text {TML}} = \mathbb {E}(k\ell _0) = \sum _{k=1}^N k\ p_k = \sum _{k=1}^N k \frac{g_k e^{-\beta k \psi _0}}{Z} = -\frac{1}{\psi _0}\frac{\partial \ln Z}{\partial \beta } = \Phi \end{aligned}$$

(17)

That is, Eq. (16) is nothing more than the statement that a higher initial quenched-disorder parameter $\beta$ will make the total microcracked length (TML) and grow faster the higher the $\beta$. And therefore, we have that in elderly people the microcracking will grow faster than in young people and that would explain loss of bone resistance with age, which is precisely what is observed in microphotographs by [63].

Moreover, the sharp decrease in the work to fracture observed in elderly people [58, 59], could be interpreted as the existence of much more extensive initial microcracking in elderly people (higher $\beta$); in fact, a significant correlation (p-value $< 0.01$) was found between $\beta$ and the factor estimating softening by microcracking w. In addition, the higher is $\beta$, the higher is the final increase in entropy $\Delta H$ (p-value $\approx 0.01$). All these facts are consistent with the initial assumptions of the model.

The proposed model was taking into account the prediction of the possible configurations of a microcrack by means of the function $g(E_k)=r^k$, being r the lattice parameter and $1\le k \le N$ the possible energy levels $E_k$ directly associated with the length of the microcrack (here assumed to be computed from the number of osteon limits crossed). The mean value of r was 1.31±0.22, within the range $1<r\le 3$ expected from the three idealized nets presented. This value is lower than those of the hexagonal lattice, which is the most tortuous net of those proposed, suggesting that the propagation of a microcrack through cortical bone is still more difficult than this idealized lattice. This could be related to the irregularity of the cortical bone microstructure, the variation of osteon diameters and the lower energy needed to propagate through certain paths that would be preferably crossed by the microcrack. Moreover, N represents the number of propagation stages that the crack follows during its progress, related to the number of osteons that have place in the net, and due to the finite thickness of cortical bone, the number of osteons crossed by the microcrack was expected to be low. The mean value of N (3.08±0.61) is consistent with the previous discussion, as the cross section of the cortical bone coupons was approximately 0.5 mm thick and 2 mm width, and the mean diameter of an osteon is 0.20 mm.

On the other hand, in its current form, this study does not provide an explanation for the observed variance in the parameters $(\mu _i, \nu _i, \rho _i, \tau _i, \varsigma _i)$ of equation (6). The relationship between these parameters and measurable densimetric magnitudes is a desirable goal for future research [64]. Due to the limited experimental design and lack of specific data, this study leaves unanswered a substantial number of potentially intriguing questions. Adding densimetric measurements to the mechanical data could have certainly clarified the reason for the observed variations in the model’s parameters.

A limitation of this study is that it focuses exclusively on understanding the effect of microcracking, so other important factors such as local bone mineral density [65,66,67], density distribution [68] and its effect on mechanical properties are left out. These are important factors, which should be taken into account in a more detailed model of the cortical bone [69,70,71]. Similarly, it would be desirable to use a generalization of this model in cases of nonmonotonic load and where the degree of microcracking could be determined by direct observation using micrographs [43, 72, 73] or other alternative techniques such as acoustic emission [40].

Finally, despite the fact that other studies have demonstrated that the ultimate stress of rupture of biological tissues is notoriously influenced by strain rate [74,75,76,77], our model provides no insight into the strain rate and its effect on the ease with which microcracking forms. In fact, in rib bone the traumatic fracture usually happens with much faster loading rates (as a nonweight bearing site, rib bone fracture mechanism and nanostructure might essentially be different from the clinically relevant site such as femoral neck). For these reasons, a viscoelastic generalization of the model, following the approach of other authors who have extensively studied such aspects [18, 78, 79], may clarify the effect of strain rate in these situations.

Appendix

1.1 Stress–strain constitutive relationships

This section presents the computational details for stress–strain relationship (10) and the entropy equation. The former is obtained by deriving the SEDF $\Psi (\beta ,\textbf{E})$ with respect to the strain tensor $\textbf{E}$ as stated in Eq. (9). Taking into account the form of $\Psi$ given by expression (2), where the partition function Z is given by (8), the explicit calculation of the derivative of SEDF with respect to the strain tensor is:

$$\begin{aligned} \textbf{S} = -\frac{1}{\beta Z}\frac{\partial Z}{\partial \textbf{E}} = -\frac{1}{\beta Z}\sum _k g_k(-\beta k)\frac{\partial \psi _0(\textbf{E})}{\partial \textbf{E}} e^{-\beta \psi _0(\textbf{E})k} = \frac{1}{Z}\frac{\partial \psi _0}{\partial \textbf{E}} \sum _k kg_k e^{-\beta \psi _0k} \end{aligned}$$

(18)

where $\partial \psi _0/\partial \textbf{E}=\bar{\textbf{S}}$ is the mechanical stress without microcracking. However, for computational purposes, it is more interesting to calculate the summation as a closed-form expression. For this purpose, we derive $\ln Z$ with respect to $\beta$ instead of $\textbf{E}$:

$$\begin{aligned} \frac{\partial \ln Z}{\partial \beta } = \frac{1}{Z}\frac{\partial Z}{\partial \beta } = \frac{1}{Z}\sum _k g_k (-k\psi _0)e^{-\beta \psi _0k}= -\frac{1}{Z}\psi _0\sum _k kg_ke^{-\beta \psi _0k} \end{aligned}$$

(19)

Note that this last relationship (19) differs from (18) in a factor $-\bar{\textbf{S}}/\psi _0$, such that:

$$\begin{aligned} \textbf{S} = -\frac{1}{\beta }\frac{\partial \ln Z}{\partial \textbf{E}} = -\frac{1}{\psi _0}\frac{\partial \ln Z}{\partial \beta }\ \bar{\textbf{S}} \end{aligned}$$

(20)

obtaining then the same result of (18) by deriving with respect to $\beta$ instead of $\textbf{E}$. This allows us to write the explicit form of $\textbf{S}$ in closed form using Z in the closed form of (8), such that:

$$\begin{aligned} \textbf{S}=-\frac{\bar{\textbf{S}}}{\psi _0}\frac{\partial }{\partial \beta } \left[ \ln \left( re^{-\beta \psi _0}\frac{1-(re^{-\beta \psi _0})^N}{1-re^{-\beta \psi _0}}\right) \right] = \bar{\textbf{S}}\left( \frac{1}{1-re^{-\beta \psi _0}}-\frac{Nr^Ne^{-\beta \psi _0 N}}{1-r^Ne^{-\beta \psi _0N}}\right) \end{aligned}$$

(21)

which is the stress–strain relationship sought, which has the general form $\textbf{S} = \bar{\textbf{S}} \Phi (\beta ,\textbf{E}),N,r)$, with $\Phi$ being the function that corrects the “uncracked stress” to include the effect of crack propagation on the mechanical behavior.

The components of the “uncracked stress” can be obtained deriving the function $\psi _0$, namely $\bar{\textbf{S}} =\partial \psi _0/\partial \textbf{E}$. The explicit derivations are:

$$\begin{aligned} \bar{S}_{xx}&= \left( (2\alpha _1+\gamma _1)I_1+3\alpha _2I_1^2+\gamma _2(I_1^2+I_2)\right) -(\gamma _1+ \gamma _2I_1)E_{xx} + \nu E_{zz}E_{yy} + \dots \nonumber \\&\quad \dots +\left( 2\kappa _1I_4+3\kappa _2I_4^2+\varsigma _2I_5\right) a_x^2+ 2\left( \varsigma _1+\varsigma _2I_4\right) (E_{xx}a_x+E_{xy}a_y)a_x \end{aligned}$$

(22a)

$$\begin{aligned} \bar{S}_{yy}&= \left( (2\alpha _1+\gamma _1)I_1+3\alpha _2I_1^2+\gamma _2(I_1^2+I_2)\right) -(\gamma _1+ \gamma _2I_1)E_{yy} + \nu E_{zz}E_{xx} + \dots \nonumber \\&\quad \dots +\left( 2\kappa _1I_4+3\kappa _2I_4^2+\varsigma _2I_5\right) a_y^2+2\left( \varsigma _1+ \varsigma _2I_4\right) (E_{xy}a_x+E_{yy}a_y)a_y=0 \end{aligned}$$

(22b)

$$\begin{aligned} \bar{S}_{zz}&= \left( (2\alpha _1+\gamma _1)I_1+3\alpha _2I_1^2+\gamma _2(I_1^2+I_2)\right) -\dots \nonumber \\&\quad \dots -(\gamma _1+ \gamma _2I_1)E_{zz} + \nu (E_{xx}E_{yy}-E_{xy}^2)=0 \end{aligned}$$

(22c)

$$\begin{aligned} \bar{S}_{xy}&= -(\gamma _1+ \gamma _2I_1)E_{xy} - \nu E_{zz}E_{xy} + \left( 2\kappa _1I_4+ 3\kappa _2I_4^2+\varsigma _2I_5\right) a_x a_y + \dots \nonumber \\&\quad \dots + \left( \varsigma _1+\varsigma _2I_4\right) (E_{xy} +(E_{xx}+E_{yy})a_x a_y) \end{aligned}$$

(22d)

These expressions are used for equations (14) and (28).

1.2 Explicit forms of strain and stress tensor for tensile and bending tests

In both types of tests, the specimen is subjected to a plane stress loading case. For the tensile tests, the Green–Lagrangian strain tensor and the second Piola–Kirchhoff stress tensor take the forms:

$$\begin{aligned} \textbf{E}(t) = \begin{bmatrix} \varepsilon _t &{} 0 &{} 0 \\[1em] 0 &{} -\bar{\nu }(\varepsilon _t)\varepsilon _t &{} 0 \\[1em] 0 &{} 0 &{} -\bar{\nu }(\varepsilon _t)\varepsilon _t \end{bmatrix} \qquad \textbf{S}(t) = \begin{bmatrix} \sigma _t &{} 0 &{} 0 \\[1em] 0 &{} 0 &{} 0 \\[1em] 0 &{} 0 &{} 0 \end{bmatrix} \end{aligned}$$

(23)

where $\bar{\nu }(\cdot )$ is a function that accounts for the Poisson’s effect (for a nonlinear material the Poisson’s ratio is no longer constant). For the uniaxial tensile test we have that the longitudinal stress is $\sigma _t =F_t/(\lambda _t A_0)$, where $F_t$ is the tensile force, $\lambda _t$ the longitudinal stretch, and $A_0$ the initial cross-sectional area. For the bending test and using curvilinear Frenet-Serret coordinates $(\xi , \eta , \zeta )$ [60], see Fig. 8, then strain and tensors take the forms:

Fig. 8

Curvilinear Frenet-Serret coordinates for the points of a rib in the bending plane

$$\begin{aligned} \textbf{E}(t) = \begin{bmatrix} \varepsilon _t &{} \gamma _t &{} 0 \\[1em] \gamma _t &{} -\bar{\nu }(\varepsilon _t)\varepsilon _t &{} 0 \\[1em] 0 &{} 0 &{} -\bar{\nu }(\varepsilon _t)\varepsilon _t \end{bmatrix} \qquad \textbf{S}(t) = \begin{bmatrix} \sigma _t &{} \tau _t &{} 0 \\[1em] \tau _t &{} 0 &{} 0 \\[1em] 0 &{} 0 &{} 0 \end{bmatrix} \end{aligned}$$

(24)

Because both $\varepsilon _t$ and $\gamma _t$ are proportional to the applied load, it turns out that we have $\gamma _t \propto \varepsilon _t$. For long bones as the human ribs, the tangential stress $\tau _t$ is much lower than the normal stress $\sigma _t$ (indeed, $\tau _t < 0.07\sigma _t$ [60]), so they play a minor role in microcracking and fracture. In addition, we have that the longitudinal strain is given by:

$$\begin{aligned} \varepsilon _t = \frac{1}{2} \left[ 1-\frac{(1-\eta \chi _0)^2}{(1-\Lambda _t\eta \chi _t)^2} \right] \end{aligned}$$

(25)

For small strain, we obtain that $\varepsilon _t \propto (\chi _t-\chi _0)$, so for bending tests the strain is initially proportional to the change of curvature. In the above formula, $\Lambda _t$ is a factor taking into account the cross-section distortion occurred during the deformation process that can be computed as [60]:

$$\begin{aligned} \Lambda _t = \frac{1-\chi \eta }{1-\chi _0\eta }\left[ 1+(1+\bar{\nu })\left( \frac{(1-\chi \eta )^2}{(1-\chi _0\eta )^2}-1\right) \right] ^{-1/2} \end{aligned}$$

(26)

1.3 Computation of the bending moments

In this section we explain the computation for the numerical value of the bending moment $\bar{M}_b^{(t)}$ and for the bending moment as a function of the model’s parameters $M_b(\chi _t;\varvec{\theta },\beta ,r,N)$, both appearing in Eq. (15).

The value of the bending moment is given by $\bar{M}_b^{(t)} = R_L a_t = F_t a_t b_t/(a_t+b_t)$ as can be seen in Fig. 9. On the other hand, the bending moment as a function of the parameter requires the computation of the following integral over the cross-section $\Omega$:

$$\begin{aligned} M_b(\chi _t,\varvec{\theta },\beta ,r,N)= \iint _\Omega \hat{\sigma }(\chi _t,\varvec{\theta },\beta ,r,N;\eta )\eta \ \text {d}\eta \text {d}\zeta \end{aligned}$$

(27)

where the function $\hat{\sigma }$ is defined as:

$$\begin{aligned} \hat{\sigma }(\chi _t,\varvec{\theta },\beta ,r,N;\eta ) = \left( \frac{1}{1-re^{-\beta \psi _0}} - \frac{Nr^Ne^{-\beta \psi _0N}}{1-r^Ne^{-\beta \psi _0N}}\right) \cdot \bar{S}_{xx}\left( \varepsilon _t(\chi _t,\eta ),\varvec{\theta }\right) \end{aligned}$$

(28)

being $\bar{S}_{xx}$ the first function in Eq. (22), and $\varepsilon _t(\chi _t,\eta )$ the strain given by Eq. (25).

Fig. 9

Calculation diagram for the bending moment, indicating the relevant reactions and distances: $a_t$ right horizontal distance, $b_t$ left horizontal distance

1.4 Values for the parameters of the model for tensile and bending tests

In this section, the mean and deviation values for the different parameters, both for tensile and bending tests, are shown, divided in age groups. Note that in both Tables 2 and 3, the value $\gamma _3 (= \rho _3 + \tau _3)$ appears, as it is not possible to distinguish the effects of the terms $\rho _3 I_4^3$ and $\tau _3 I_4^2 I_5$ of equation (6) in tensile and bending tests. To differentiate between the two parameters, shear or torsion tests would have been necessary, which was beyond the scope of this study.

Table 2 Parameter values for the coupon specimens under uniaxial tensile tests fitted using the model proposed, in age groups

Table 3 Parameter values for the bending specimens fitted using the model proposed, in age groups

In general, in Tables 2 and 3, both the sign and order of magnitude of the parameters are consistent, with slightly larger errors for bending tests. In particular, the parameters associated with second-order terms in the strain $(\mu _2, \nu _2, \rho _2, \tau _2)$ as well as the crack parameters $(\beta , N, r)$ obtained for tensile test are very similar to those obtained for bending tests for all age groups. On the other hand, it is also observed that the third-order parameters ($\square _3$) have larger error estimates, and it is precisely in two of these parameters, $\gamma _3$ and $\nu _3$, where sign differences are observed for certain age groups. We believe that these discrepancies are the result of a greater estimation error in the cubic terms of strain Eq. (6).