Analysis of Damage Models for Cortical Bone

Abstract

Bone tissue is a material with a complex structure and mechanical properties. Diseases or even normal repetitive loads may cause microfractures to appear in the bone structure, leading to a deterioration of its properties. A better understanding of this phenomenon will lead to better predictions of bone fracture or bone-implant performance. In this work, the model proposed by Frémond and Nedjar in 1996 (initially for concrete structures) is numerically analyzed and compared against a bone specific mechanical model proposed by García et al. in 2009. The objective is to evaluate both models implemented with a finite element method. This will allow us to determine if the modified Frémond–Nedjar model is adequate for this purpose. We show that, in one dimension, both models show similar results, reproducing the qualitative behaviour of bone subjected to typical engineering tests. In particular, the Frémond–Nedjar model with the introduced modifications shows good agreement with experimental data. Finally, some two-dimensional results are also provided for the Frémond–Nedjar model to show its behaviour in the simulation of a real tensile test.

1. Introduction

Damage models arise in order to describe how mechanical properties of materials degrade over time. This degradation can be caused both because of the loading it is subjected to, and due to external causes (such as crack formation due to thermal shock or chemical attack). These models have been deeply studied for structures, usually concrete structures, where the progressive wear of the materials can be critical to its integrity. In this field, works about damage exist since the decade of 1980 [1]. A later work from Frémond and Nedjar in 1996 [2] became a reference for new concrete damage models based on the continuum elastoplastic damage approach [3]. In this approach, the principle of virtual power is modified, including the damage in the term of the power of the internal forces. Also, damage is represented by a scalar field.

The study of engineering concepts in the field of biology and medicine is more recent, but it is growing fast. In the particular case of bone tissue, biological aspects such as bone remodelling were studied also since the decade of 1980 [4]. The model proposed by Weinans, Huiskes and Grootenboer in 1992 [5] started the possibility of simulating this effect numerically, and led to a large number of contributions in this field. The effect of damage described previously can be seen in bones too. In bone tissue, both loading and external causes (now related to illness such as osteoporosis) produce again a growing deterioration of its elastic properties [6]. The first studies of damage were focused on cumulative damage caused by cyclic loading [7,8]. This approach allows for fatigue estimations of the number of cycles that a probe can withstand, but it is not effective for its implementation in numerical simulations. In 1999, Fondrk, Bahniuk and Davy proposed a model that reproduced the tensile behaviour of cortical bone [9]; however, it is limited to one-dimensional simulations, obtaining bending results by applying beam theory assumptions. A recent work from and García et al. in 2009 [10] introduced rheological bone specific models that reproduced uniaxial and cyclic tests and could be extended to three dimensional simulations. Other works about bone damage can be seen in [11,12,13,14].

In this work, the Frémond–Nedjar model is compared against the García et al. model in order to assess its capabilities to reproduce bone tissue behaviour. Indeed, one of the main novelties of this work is its application in the simulation of the bone damage process. The selection and the interest of analysing the Frémond-Nejdar model relies on the fact that its formulated as a set of partial differential equations (in particular a subdifferential inclusion) that couples the damage evolution with the well known linear elasticity model. This allows for its numerical analysis and its implementation in a finite element simulation. Also, it is not restricted to one dimensional simulations, like the previously mentioned Fondrk, Bahniuk and Davy model.

Furthermore, since the formulation is similar to the Weinans–Huiskes–Grootenboer model of bone remodelling, it would allow for a direct coupling of these models, a future objective of the authors. Although an existing work presents a coupling between damage and a remodelling model [15], it is based on a simple remodelling rule and the numerical analysis is not performed.

The paper is structured as follows. In Section 2, the formulations of the studied models are presented, then, in Section 3, the implementation and results obtained are shown and discussed, followed by some conclusions in Section 4.

2. Damage Models

In this paper, as mentioned before, two damage models for numerical simulations are analyzed: the Frémond–Nedjar model, first developed for concrete structures, and the model proposed by García et al. (bone specific). The particular formulation of each model is described in the following subsections. Special emphasis is placed on the Frémond–Nedjar model, for which the numerical analysis of the algorithm proposed for its resolution is shown.

2.1. Frémond–Nedjar Model

The damage model proposed by Frémond and Nedjar in 1996 considers damage as an unknown of the problem ($\beta$) that varies between 1 (undamaged material) and 0 (completely damaged). This counter-intuitive definition accounts for the fact that the variable multiplies the mechanical properties of the material (the elastic modulus in the one-dimensional case) in such a way that, in a damaged material, they will be affected by this variable.

As mentioned before, this model was presented for concrete structures, so some modifications to include bone behaviour were necessary. Since bone tissue is a living material that can heal over time, the assumption that this variable cannot recover, $\dot{\beta }<0$, is no longer true and it is removed from the initial formulation.

The mathematical formulation of this model with the modifications to account for bone properties is defined in what follows.

In [2] the damage source function $\varphi$ was defined by

$$\varphi (\mathit{\epsilon }\left(\mathit{u}\right),\beta )={\lambda }_d\left(\frac{1-\beta }{\beta }\right)-{\lambda }_u{\left|\mathit{\epsilon }\left(\mathit{u}\right)\right|}^2+{\lambda }_w,$$

where ${\lambda }_d,{\lambda }_u$ and ${\lambda }_w$ are constitutive parameters. The second term becomes unmanageable mathematically when strains are very large, but then the whole model becomes inadequate, so we truncate the term as follows. Given $q^{\ast }>0$, a sufficiently large strain energy truncation constant, let ${\mathsf{\Psi }}_{q^{\ast }}$ be given by

$${\mathsf{\Psi }}_{q^{\ast }}\left(\mathit{\tau }\right)=\left\{\begin{array}{c}{\left|\mathit{\tau }\right|}^2\mathrm{if}{\left|\mathit{\tau }\right|}^2\le q^{\ast },\hfill \\ q^{\ast }\mathrm{otherwise},\hfill \end{array}\right.$$

where ${\left|\mathit{\tau }\right|}^2={\tau }_{ij}{\tau }_{ij}$, and here and below, $i,j=1,\dots ,d$ (d is the dimension of the problem), and a repeated index indicates summation. Therefore, we use the truncated damage source function:

$$\varphi (\mathit{\epsilon }\left(\mathit{u}\right),\beta )={\lambda }_d\left(\frac{1-\beta }{\beta }\right)-{\lambda }_u{\mathsf{\Psi }}_{q^{\ast }}\left(\mathit{\epsilon }\left(\mathit{u}\right)\right)+{\lambda }_w.$$

(1)

Finally, the mechanical problem is written as follows.

Problem P. Find the displacement field $\mathit{u}={\left(u_i\right)}_{i=1}^d:\overline{\mathsf{\Omega }}\times [0,T]\to {\mathbb{R}}^d$ and the damage field $\beta :\overline{\mathsf{\Omega }}\times [0,T]\to \mathbb{R}$ such that

$$\begin{array}{c}\mathit{\sigma }=2\beta \mu \mathit{\epsilon }\left(\mathit{u}\right)+\beta \lambda \mathrm{Div}\left(\mathit{u}\right)\mathcal{I}\mathrm{in}\mathsf{\Omega }\times (0,T),\hfill \end{array}$$

(2)

$$\begin{array}{c}-\mathrm{Div}\mathit{\sigma }={\mathit{f}}_0\mathrm{in}\mathsf{\Omega }\times (0,T),\hfill \end{array}$$

(3)

$$\begin{array}{c}\dot{\beta }-{\kappa }^{\ast }\Delta \beta +\partial I_{[{\beta }_{\ast },1]}\left(\beta \right)\ni \varphi (\mathit{\epsilon }\left(\mathit{u}\right),\beta )\mathrm{in}\mathsf{\Omega }\times (0,T),\hfill \end{array}$$

(4)

$$\begin{array}{c}\mathit{u}=\mathbf{0}\mathrm{on}{\mathsf{\Gamma }}_D\times (0,T),\hfill \end{array}$$

(5)

$$\begin{array}{c}\frac{\partial \beta }{\partial \mathit{\nu }}=0\mathrm{on}\mathsf{\Gamma }\times (0,T),\hfill \end{array}$$

(6)

$$\begin{array}{c}\mathit{\sigma }\mathit{\nu }=\mathit{g}\mathrm{on}{\mathsf{\Gamma }}_N\times (0,T),\hfill \end{array}$$

(7)

$$\begin{array}{c}\beta (\mathit{x},0)={\beta }_0\left(\mathit{x}\right)\mathrm{in}\mathsf{\Omega },\hfill \end{array}$$

(8)

where $\mathsf{\Omega }$ represents the bone, whose boundary $\mathsf{\Gamma }$ is assumed to be decomposed into two parts ${\mathsf{\Gamma }}_D$ and ${\mathsf{\Gamma }}_N$ such that $\mathsf{\Gamma }={\mathsf{\Gamma }}_D\cup {\mathsf{\Gamma }}_N$, with $meas\left({\mathsf{\Gamma }}_D\right)>0$, $[0,T]$, for $T>0$, denotes the time interval of interest, and let $\mathrm{Div}$ be the divergence of tensor-valued functions. $\lambda$ and $\mu$ are the classical Lame’s coefficients and ${\kappa }^{\ast }>0$ is a damage diffusion coefficient. Moreover, ${\beta }_0$ is an initial condition for the damage field, ${\mathit{f}}_0$ and g represent volume and traction forces, respectively, and we include the subdifferential of the indicator function $I_{[{\beta }_{\ast },1]}$ into (4) to ensure that the damage function belongs to the interval $[{\beta }_{\ast },1]$. Here, ${\beta }_{\ast }>0$ is a positive constant, assumed small, and it is introduced for mathematical reasons. In any case, when damage becomes zero, the material is dense with microcracks and modelling it as elastic ceases to make sense (see [16] for details). Finally, we note that damage function $\varphi$ is given in (1), where constants ${\lambda }_d$, ${\lambda }_u$ and ${\lambda }_w$ are constitutive parameters.

Now, in order to obtain the variational formulation of Problem P, let $Y=L^2\left(\mathsf{\Omega }\right)$, $H={\left[L^2\left(\mathsf{\Omega }\right)\right]}^d$ and $Q={\left[L^2\left(\mathsf{\Omega }\right)\right]}^{d\times d}$ and denote by ${(·,·)}_Y$, ${(·,·)}_H$ and ${(·,·)}_Q$ the respective scalar products in these spaces. Moreover, let us define the variational space V as follows,

$$\begin{array}{c}V=\{\mathit{v}\in {\left[H^1\left(\mathsf{\Omega }\right)\right]}^d;\mathit{v}=\mathbf{0}\mathrm{on}{\mathsf{\Gamma }}_D\},\hfill \end{array}$$

with the scalar product ${(·,·)}_V$, and norm ${\parallel ·\parallel }_V$.

Finally, let us define the convex set of admissible damage functions,

$$\begin{array}{c}K=\{\zeta \in H^1\left(\mathsf{\Omega }\right);{\beta }_{\ast }\le \zeta \le 1\mathrm{a}.\mathrm{e}. \mathrm{in} \mathsf{\Omega }\}.\hfill \end{array}$$

By using Green’s formula and boundary conditions (5)–(7), we write the variational formulation of problem P.

Problem VP. Find the displacement field $\mathit{u}:[0,T]\to V$ and the damage field $\beta :[0,T]\to K$ such that $\beta \left(0\right)={\beta }_0$ and, for a.e. $t\in (0,T)$, and for all $\mathit{v}\in V$ and $\zeta \in K$,

$$\begin{array}{c}c(\beta \left(t\right);\mathit{u}\left(t\right),\mathit{v})={(\mathit{f}\left(t\right),\mathit{v})}_V,\hfill \end{array}$$

(9)

$$\begin{array}{c}{\displaystyle {(\dot{\beta }\left(t\right),\zeta -\beta \left(t\right))}_Y+{\kappa }^{\ast }{(\nabla \beta \left(t\right),\nabla (\zeta -\beta \left(t\right)))}_H\ge {\left(\varphi \left(\mathit{\epsilon }\right(\mathit{u}\left(t\right)),\beta (t\left)\right),\zeta -\beta \left(t\right)\right)}_Y,}\hfill \end{array}$$

(10)

where the bilinear functional $c:K\times V\times V\to \mathbb{R}$ is defined as follows, for all $\mathit{v},\mathit{w}\in V$ and $\zeta \in K$,

$$c(\zeta ;\mathit{v},\mathit{w})={\int }_{\mathsf{\Omega }}\zeta \left(2\mu \mathit{\epsilon }\left(\mathit{v}\right):\mathit{\epsilon }\left(\mathit{w}\right)+\lambda \mathrm{div}\mathit{v}\mathrm{div}\mathit{w}\right)d\mathit{x}.$$

The operator div represents now the divergence of vector-valued functions, the element $\mathit{f}\left(t\right)\in V^{\prime }$ (as usual, $V^{\prime }$ denotes the dual space of V) is obtained from Riesz’ theorem as

$${(\mathit{f}\left(t\right),\mathit{w})}_V={({\mathit{f}}_0\left(t\right),\mathit{w})}_H+{(\mathit{g}\left(t\right),\mathit{w})}_{{\left[L^2\left({\mathsf{\Gamma }}_N\right)\right]}^d}\forall \mathit{w}\in V,$$

and we recall that damage function $\varphi$ is defined in (1).

We now consider a fully discrete approximation of problem $VP$. To solve it numerically, a finite element scheme is used. This is done in two steps. First, we assume that the domain $\overline{\mathsf{\Omega }}$ is polyhedral and we denote by ${\mathcal{T}}^h$ a regular triangulation in the sense of [17]. Thus, we construct the finite dimensional spaces $V^h\subset V$ and $E^h\subset H^1\left(\mathsf{\Omega }\right)$ given by

$$\begin{array}{c}{V}^h=\{{\mathit{v}}^h\in {\left[C\left(\overline{\mathsf{\Omega }}\right)\right]}^d;{\mathit{v}}_{|Tr}^h\in {\left[P_1\left(Tr\right)\right]}^d\forall Tr\in {\mathcal{T}}^h,{\mathit{v}}^h=\mathbf{0}\mathrm{on}{\mathsf{\Gamma }}_D\},\hfill \\ E^h=\{{\xi }^h\in C\left(\overline{\mathsf{\Omega }}\right);{\xi }_{|Tr}^h\in P_1\left(Tr\right)\forall Tr\in {\mathcal{T}}^h\},\hfill \end{array}$$

where $P_1\left(Tr\right)$ represents the space of polynomials of degree less or equal to one in the element $Tr$, i.e., the finite element spaces $V^h$ and $E^h$ are composed of continuous and piecewise affine functions. Here, $h>0$ denotes the spatial discretization parameter. Moreover, let $K^h=K\cap E^h$ and assume that the discrete initial condition, denoted by ${\beta }_0^h$, is given by

$${\beta }_0^h={\mathcal{P}}^h{\beta }_0,$$

where ${\mathcal{P}}^h$ is the classical finite element interpolation operator over $E^h$ (see, e.g., [17]).

Secondly, we consider a partition of the time interval $[0,T]$, denoted by $0=t_0<t_1<\cdots <t_N=T$. In this case, we use a uniform partition with step size $k=T/N$ and nodes $t_n=nk$ for $n=0,1,\dots ,N$.

Therefore, using a combination of both implicit and explicit Euler schemes, the fully discrete approximations are considered as follows.

Problem VP${}^{hk}$. Find the discrete displacement field ${\mathit{u}}^{hk}={\left\{{\mathit{u}}_n^{hk}\right\}}_{n=0}^N\subset V^h$ and the discrete damage field ${\beta }^{hk}={\left\{{\beta }_n^{hk}\right\}}_{n=0}^N\subset K^h$ such that ${\beta }_0^{hk}={\beta }_0^h$ and, for all ${\mathit{v}}^h\in V^h$, ${\zeta }^h\in K^h$,

$$\begin{array}{c}c({\beta }_n^{hk};{\mathit{u}}_n^{hk},{\mathit{v}}^h)={({\mathit{f}}_n,{\mathit{v}}^h)}_{V}n=0,1,\dots ,N,\hfill \\ {\displaystyle {({\beta }_n^{hk},{\zeta }^h-{\beta }_n^{hk})}_Y+k{\kappa }^{\ast }{(\nabla {\beta }_n^{hk},\nabla ({\zeta }^h-{\beta }_n^{hk}))}_H\ge {({\beta }_{n-1}^{hk},{\zeta }^h-{\beta }_n^{hk})}_Y}\hfill \end{array}$$

(11)

$$\begin{array}{c}{\displaystyle +k{\left(\varphi (\mathit{\epsilon }\left({\mathit{u}}_{n-1}^{hk}\right),{\beta }_{n-1}^{hk}),{\zeta }^h-{\beta }_n^{hk}\right)}_{Y}n=1,\dots ,N.}\hfill \end{array}$$

(12)

We remark that in problem VP${}^{hk}$ the “initial condition” ${\mathit{u}}_0^{hk}$ for the displacement field must be calculated because it is not previously given. Thus, we take it as the solution to the corresponding discrete problem:

$$c({\beta }_0^h;{\mathit{u}}_0^{hk},{\mathit{v}}^h)={({\mathit{f}}_0,{\mathit{v}}^h)}_V\forall {\mathit{v}}^h\in V^h.$$

We have the following theorem which states the linear convergence of these approximations under suitable additional regularity conditions.

Theorem 1.

Assume that problem VP has a unique solution $(\mathit{u},\beta )$ with the following regularity:

$$\begin{array}{c}\mathit{u}\in \mathcal{C}([0,T];{\left[H^2\left(\mathsf{\Omega }\right)\right]}^d)\cap {\mathcal{C}}^1([0,T];V),\hfill \\ \beta \in H^2(0,T;Y)\cap H^1(0,T;H^1\left(\mathsf{\Omega }\right))\cap \mathcal{C}([0,T];H^2\left(\mathsf{\Omega }\right)),\hfill \end{array}$$

and let $({\mathit{u}}^{hk},{\beta }^{hk})$ be the solution to Problem VP${}^{hk}$. Then, it follows that the numerical approximation is linearly convergent; that is, there exists a positive constant C, independent of the discretization parameters h and k, such that

$$\underset{0\le n\le N}{max}\left\{\parallel {\mathit{u}}_n-{\mathit{u}}_n^{hk}{\parallel }_V+{\parallel {\beta }_n-{\beta }_n^{hk}\parallel }_Y\right\}\le C(h+k).$$

Proof. 

First, proceeding as in [16], we have the following estimates for the damage field, for all ${\zeta }^h={\{{\zeta }_j^h\}}_{j=0}^n\subset K^h$,

$$\begin{array}{c}{\displaystyle \parallel {\beta }_n-{\beta }_n^{hk}{\parallel }_Y^2+k\sum _{j=1}^n\parallel \nabla ({\beta }_n-{\beta }_n^{hk}){\parallel }_Y^2\le C(\parallel {\beta }_0-{\beta }_0^h{\parallel }_Y^2+\parallel {\beta }_1-{\zeta }_1^h{\parallel }_Y^2+{\parallel {\beta }_n-{\zeta }_n^h\parallel }_Y^2}\hfill \\ {\displaystyle +k\sum _{j=1}^n\parallel {\beta }_{j-1}-{\beta }_{j-1}^{hk}{\parallel }_Y^2+k\sum _{j=1}^n\parallel {\mathit{u}}_j-{\mathit{u}}_j^{hk}{\parallel }_V^2+k\sum _{j=1}^n\parallel {\beta }_j-{\zeta }_j^h{\parallel }_{H^1\left(\mathsf{\Omega }\right)}^2+k^2+k\sum _{j=1}^n{\parallel \delta {\beta }_j-{\dot{\beta }}_j\parallel }_Y^2}\hfill \\ {\displaystyle +k\sum _{j=1}^n\parallel \varphi (\mathit{\epsilon }\left({\mathit{u}}_j\right),{\beta }_j)-\delta {\beta }_j+{\kappa }^{\ast }\Delta {\beta }_j{\parallel }_Y\parallel {\beta }_j-{\zeta }_j^h{\parallel }_Y+k^{-1}\sum _{j=1}^{n-1}{\parallel {\beta }_{j+1}-{\zeta }_{j+1}^h-({\beta }_j-{\zeta }_j^j)\parallel }_Y^2),}\hfill \end{array}$$

where we used the notation $\delta {\beta }_j=({\beta }_j-{\beta }_{j-1})/k$ and the regularity of the continuous solution $\beta$.

Now, we obtain the estimates for the displacement fields. Subtracting the variational Equation (9) at time $t_n$ for $\mathit{v}={\mathit{v}}^h\in V^h\subset V$ and the discrete variational Equation (11) we have

$$c({\beta }_n;{\mathit{u}}_n,{\mathit{v}}^h)-c({\beta }_n^{hk};{\mathit{u}}_n^{hk},{\mathit{v}}^h)\forall {\mathit{v}}^h\in V^h,$$

so it follows that, for all ${\mathit{v}}^h\in V^h,$

$${\displaystyle c({\beta }_n;{\mathit{u}}_n,{\mathit{u}}_n-{\mathit{u}}_n^{hk})-c({\beta }_n^{hk};{\mathit{u}}_n^{hk},{\mathit{u}}_n-{\mathit{u}}_n^{hk})=c({\beta }_n;{\mathit{u}}_n,{\mathit{u}}_n-{\mathit{v}}^h)-c({\beta }_n^{hk};{\mathit{u}}_n^{hk},{\mathit{u}}_n-{\mathit{v}}^h).}$$

Taking into account that

$$\begin{array}{c}{\displaystyle c({\beta }_n;{\mathit{u}}_n,{\mathit{u}}_n-{\mathit{u}}_n^{hk})-c({\beta }_n^{hk};{\mathit{u}}_n^{hk},{\mathit{u}}_n-{\mathit{u}}_n^{hk})=c({\beta }_n^{hk};{\mathit{u}}_n-{\mathit{u}}_n^{hk},{\mathit{u}}_n-{\mathit{u}}_n^{hk})+c({\beta }_n-{\beta }_n^{hk};{\mathit{u}}_n,{\mathit{u}}_n-{\mathit{u}}_n^{hk}),}\hfill \\ c({\beta }_n^{hk};{\mathit{u}}_n-{\mathit{u}}_n^{hk},{\mathit{u}}_n-{\mathit{u}}_n^{hk})\ge C{\parallel {\mathit{u}}_n-{\mathit{u}}_n^{hk}\parallel }_V^2,\hfill \end{array}$$

using the fact that ${\beta }_n^{hk}\in K^h$ (and so, ${\beta }_n^{hk}\ge {\beta }_{\ast }$) and the regularity $\mathit{u}\in \mathcal{C}([0,T];{\left[H^2\left(\mathsf{\Omega }\right)\right]}^d)$ we have

$$\parallel {\mathit{u}}_n-{\mathit{u}}_n^{hk}{\parallel }_V^2\le C\left(\parallel {\beta }_n-{\beta }_n^{hk}{\parallel }_Y^2+{\parallel {\mathit{u}}_n-{\mathit{v}}^h\parallel }_V^2\right)\forall {\mathit{v}}^h\in V^h.$$

Therefore, combining the previous estimates of both damage and displacement fields we conclude that

$$\begin{array}{c}{\displaystyle \parallel {\beta }_n-{\beta }_n^{hk}{\parallel }_Y^2+\parallel {\mathit{u}}_n-{\mathit{u}}_n^{hk}{\parallel }_V^2\le C(\parallel {\beta }_0-{\beta }_0^h{\parallel }_Y^2+\parallel {\beta }_1-{\zeta }_1^h{\parallel }_Y^2+\parallel {\beta }_n-{\zeta }_n^h{\parallel }_Y^2+k\sum _{j=1}^n{\parallel {\beta }_{j-1}-{\beta }_{j-1}^{hk}\parallel }_Y^2}\hfill \\ {\displaystyle +k\sum _{j=1}^n\parallel {\mathit{u}}_j-{\mathit{u}}_j^{hk}{\parallel }_V^2+k\sum _{j=1}^n\parallel {\beta }_j-{\zeta }_j^h{\parallel }_{H^1\left(\mathsf{\Omega }\right)}^2+k^2+k\sum _{j=1}^n\parallel \delta {\beta }_j-{\dot{\beta }}_j{\parallel }_Y^2+{\parallel {\mathit{u}}_n-{\mathit{v}}^h\parallel }_V^2}\hfill \\ {\displaystyle +k\sum _{j=1}^n\parallel \varphi (\mathit{\epsilon }\left({\mathit{u}}_j\right),{\beta }_j)-\delta {\beta }_j+{\kappa }^{\ast }\Delta {\beta }_j{\parallel }_Y\parallel {\beta }_j-{\zeta }_j^h{\parallel }_Y+k^{-1}\sum _{j=1}^{n-1}{\parallel {\beta }_{j+1}-{\zeta }_{j+1}^h-({\beta }_j-{\zeta }_j^j)\parallel }_Y^2).}\hfill \end{array}$$

Finally, keeping in mind that

$$k^{-1}\sum _{j=1}^{n-1}\parallel {\beta }_{j+1}-{\zeta }_{j+1}^h-({\beta }_j-{\zeta }_j^j){\parallel }_Y^2\le C{h}^2{\parallel \beta \parallel }_{H^1(0,T;H^1\left(\mathsf{\Omega }\right))}^2,$$

using the well-known approximation properties by finite elements (see [17]) and a discrete version of Gronwall’s inequality, we obtain the desired linear convergence.  □

The fully discrete scheme provided in problem VP${}^{hk}$ has been solved using a penalty-duality algorithm, related to Uzawa’s algorithm, for the numerical resolution of the discrete variational inequality. The discrete displacements have been obtained solving the discrete linear variational equation with the classical conjugate gradient method. We note that a similar scheme has also been employed for the numerical approximation of dynamic contact problems (see, for example, [16]). Moreover, the resulting algorithm has been implemented within the well-known code MATLAB in a 3.3 GHz PC (with 16 Gb of RAM memory), and a typical 1D run with parameters $h=k={10}^2$ took about $1.54$ s of CPU time.

In order to show the numerical convergence of the algorithm we solve Problem P with the following parameters:

$$\begin{array}{c}T=2,\mathsf{\Omega }=(0,1),{\mathsf{\Gamma }}_D=\left\{0\right\},{\mathsf{\Gamma }}_N=\left\{1\right\},\mathit{g}=\mathbf{0},\hfill \\ \lambda =0,\mu =1.348\times {10}^{10},{\kappa }^{\ast }=0.5,{\beta }_{\ast }=0.01,\hfill \\ q^{\ast }={10}^5,{\lambda }_d=17,{\lambda }_u=4.2\times {10}^5,{\lambda }_w=0,\hfill \end{array}$$

with the initial condition:

$${\beta }_0\left(x\right)=1\mathrm{for}\mathrm{all}x\in (0,1),$$

and the volumetric force:

$${\mathit{f}}_0(x,t)=\left\{\begin{array}{cc}90\times {10}^6·t\hfill & \hfill \mathrm{if}t<1,\\ 90\times {10}^6\hfill & \hfill \mathrm{if}t\ge 1.\end{array}\right.$$

The numerical errors are given by

$$E^{hk}=\underset{0\le n\le N}{max}\left\{\parallel {\mathit{u}}_n-{\mathit{u}}_n^{hk}{\parallel }_V+{\parallel {\beta }_n-{\beta }_n^{hk}\parallel }_Y\right\},$$

considering as “exact solution” $({\mathit{u}}_n,{\beta }_n)$ the one obtained for $h=2^{-12}$ and $k={10}^{-5}$. The errors (multiplied by 100), obtained for different discretizations, are shown in

Table 1. Numerical errors ($\times 100$) for some discretization parameters.

$\mathit{h}\downarrow \mathit{k}\to$	0.1	0.01	0.001	0.0001
$2^{-2}$	13.5490	13.5613	13.561268	13.561268
$2^{-3}$	6.6596	6.6679	6.667863	6.667863
$2^{-4}$	3.3080	3.3054	3.305386	3.305386
$2^{-5}$	1.6641	1.6457	1.645661	1.645661
$2^{-6}$	0.8460	0.8210	0.821036	0.821036
$2^{-7}$	0.4370	0.4099	0.409942	0.409942
$2^{-8}$	0.2328	0.2046	0.204565	0.204565
$2^{-9}$	0.1315	0.1017	0.101652	0.101653
$2^{-10}$	0.0822	0.0496	0.049597	0.049596

and depicted in Figure 1 against $h+k$. As shown, the linear convergence of the algorithm stated in Theorem 1 seems to be achieved.

2.2. García et al. Model

The other model examined is the one proposed by García et al. [10,18]. It consists of a rheological model developed specifically for cortical bone, to reproduce the macroscopical phenomenon observed in this tissue. The model is written using several internal variables for damage and plastic strain, as well as laws to describe the evolution of these internal variables.

In [10], several models are presented. For the present study the rate independent model is chosen, since it allows for a more immediate comparison. The rheological model is composed of an elastic spring in series with the damageable part, which consists of a secondary spring, whose elasticity varies with damage, and a friction element that determines the plasticity threshold.

From this rheological model a free energy potential is obtained, which is convex and nonsmooth:

$$\mathsf{\Psi }(\epsilon ,{\epsilon }^p,D)=\left\{\begin{array}{c}{\displaystyle \frac{1}{2}{E}_0{(\epsilon -{\epsilon }^p)}^2+\frac{1}{2}{E}_0\frac{1-D}{D}{{\epsilon }^p}^2+I_{[0,1]}\left(D\right)\mathrm{if}D>0,}\hfill \\ {\displaystyle \frac{1}{2}{E}_0{\epsilon }^2+I_{\left\{0\right\}}\left({\epsilon }^p\right)\mathrm{if}D=0,}\hfill \end{array}\right.$$

(13)

and the state laws of the material can be derived from this potential. The details of this derivation can be found in [10], and we omit them for the sake of clarity in the presentation.

Regarding the García et al. model, the algorithm they developed to solve their model was used. Again, we refer to [10] for more details about its implementation. It is based on the combination of the classical finite element method with a Newton integration scheme and a projection operator. The latter one is used to satisfy the criterion defined for the internal variables.

3. Numerical Results and Discussion

3.1. Comparison in a One-Dimensional Problem

To compare the performance of the models presented in the previous section, a one-dimensional version of both was implemented. To test the models, a typical tensile test was reproduced, since there is experimental data in the literature to evaluate the performance [18]. The one-dimensional model represents a bone fixed on one end and with an increasing displacement imposed on the other end. Stress was computed in postprocessing once the deformation was obtained.

First, we solve the Frémond–Nedjar model with the following data:

$$\begin{array}{c}T=2,\mathsf{\Omega }=(0,1),{\mathsf{\Gamma }}_D=\{0,1\},{\mathsf{\Gamma }}_N=\varnothing ,{\mathit{f}}_0=\mathbf{0},\hfill \\ \mathit{g}=\mathbf{0},\lambda =0,\mu =1.348\times {10}^{10},{\kappa }^{\ast }=0.5,{\beta }_{\ast }=0.01,\hfill \\ q^{\ast }={10}^5,{\lambda }_d=17,{\lambda }_u=4.2\times {10}^5,{\lambda }_w=0,\hfill \end{array}$$

and the initial condition:

$${\beta }_0\left(x\right)=1\mathrm{for}\mathrm{all}x\in (0,1).$$

That is, at initial time we assumed that the bone was completely healed. Moreover, since no mechanical forces were applied, the deformation is produced defining an imposed displacement at the right corner as,

$$u(1,t)=5\times {10}^{-3}·t\mathrm{for}\mathrm{all}t\in [0,2].$$

We note that the modifications needed to include this case into problem P were really straightforward, so the analysis performed in the previous section could be extended easily.

Using the discretization parameters $k=0.01$ and $h=0.05$ in Figure 2 (up) the stress–strain curve obtained with the Frémond–Nedjar model is shown. As it can be seen in the evolution of the damage variable (down), the mechanical properties of the bone degrade as the strain increases, leading to the “plastic” region which corresponds with the damage regime of the stress–strain curve.

Secondly, the García et al. model was solved, so in Figure 3 the shape of the stress–strain curve seen before is reproduced again with the rheological model.

Finally, both models were compared with experimental data obtained from [18]. Figure 4 shows the stress–strain curves obtained from both models against the experimental data (black dots). Both models show good agreement with the experimental results, with small differences in the linear part of the curve and the beginning of the damage range. The Frémond–Nedjar model shows better agreement in the linear part, but the curvature of the damage range fits worse than the rheological model. However, these small differences were not significant, since mechanical properties vary greatly from bone to bone, making very accurate fittings not useful.

In any case, both models were capable of reproducing the qualitative characteristics of bone mechanics. In particular, the modified Frémond–Nedjar model shows its capability to be used to model bone tissue in spite of it being proposed to model concrete structures.

3.2. A Two-Dimensional Problem Solved Using Frémond–Nedjar Model

Now, we consider a two-dimensional case to simulate a more realistic setting. Thus, the bone occupied a two-dimensional domain $\mathsf{\Omega }$, which was assumed fixed on its left part ${\mathsf{\Gamma }}_D$ (the whole part clamped on the vertical direction and its lower point also in the horizontal one), and subjected to the action of a surface force on its right part ${\mathsf{\Gamma }}_N$ (see Figure 5).

The following data have been used in this example:

$$\begin{array}{c}T=5,{\mathit{f}}_0=\mathbf{0},\mathit{g}=(5,0),\lambda =21.1\times {10}^{10},\mu =24.61\times {10}^{10},{\kappa }^{\ast }=1,{\beta }_{\ast }=0.01,\hfill \\ q^{\ast }={10}^5,{\lambda }_d=0.01,{\lambda }_u=15,{\lambda }_w=-0.0001,\hfill \end{array}$$

and the initial condition:

$${\beta }_0(x,y)=1\mathrm{for}\mathrm{all}(x,y)\in \mathsf{\Omega }.$$

That is, again at initial time we assumed that the bone was completely healed.

Using the time discretization parameter $k=0.1$ and the finite element mesh shown in Figure 5, the damage field at final time is plotted over the deformed configuration in Figure 6. As expected, the most damaged areas concentrated in the middle of the bone due to the applied force and the clamping conditions.

Finally, in Figure 7 the stress field is shown at final time over the deformed configuration. We note that now the highest stressed areas concentrated in the middle part of the bone. The concentration of the damage and stress in the middle of the sample agree with the observed fracture of the real probes in uniaxial tests by on middle section.

4. Conclusions

In this work, two ways of modelling damage in cortical bone were studied. The first one was a modification of a well-known damage model for concrete structures, never used for bone tissue. The second one was developed specifically for bone damage and based on a rheological model. One-dimensional simulations were performed to compare both models, reproducing a classical tensile test. Both models show similar solutions and a good agreement with experimental data. Moreover, a two-dimensional example was also considered using the Frémond–Nedjar model to show the behaviour of its solution in a more real situation. The advantage of the use of the Frémond–Nedjar model relies in its formulation, which makes it easier to couple with bone remodelling models, and it allows for a formal numerical analysis. The number of parameters required for this model is also reduced with respect to the García et al. model. These results open the possibility of using this model in bone tissue simulations.