# A gradient-enhanced bone remodelling approach to avoid the checkerboard phenomenon

Fynn Bensel, Marlis Reiber, Elise Foulatier, Philipp Junker, Udo Nackenhorst

► To cite this version:

Fynn Bensel, Marlis Reiber, Elise Foulatier, Philipp Junker, Udo Nackenhorst. A gradient-enhanced bone remodelling approach to avoid the checkerboard phenomenon. Computational Mechanics, 2023, 73, pp.1335 - 1349. 10.1007/s00466-023-02413-9. hal-05368002

HAL Id: hal-05368002

https://hal.science/hal-05368002v1

Submitted on 17 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A gradient-enhanced bone remodelling approach to avoid the checkerboard phenomenon

Fynn Bensel<fnref n="1" />,<fnref n="3" /> ID · Marlis Reiber<fnref n="1" />,<fnref n="4" /> ID · Elise Foulatier<fnref n="5" /> · Philipp Junker<fnref n="2" />,<fnref n="3" />,<fnref n="4" /> ID · Udo Nackenhorst<fnref n="1" />,<fnref n="3" />,<fnref n="4" /> ID

Received: 26 April 2023 / Accepted: 13 October 2023 / Published online: 28 November 2023

© The Author(s) 2023

## Abstract

Numerical simulation of bone remodelling enables the investigation of short- and long-term stability of bone implants and thus can be an essential tool for surgical planning. The first development of related mathematical models dates back to the early 90's, and these models have been continuously refined since then. However, one issue which has been under discussion since those early days concerns a numerical instability known as checkerboarding. A literature review of recent approaches guided us to adopt a technique established in damage mechanics and topology optimisation, where similar mesh dependencies and instabilities occur. In our investigations, the so-called gradient enhancement is used to regularise the internal variable field, representing the evolution of the bone mass density. For this, a well-established mathematical model for load-adaptive bone remodelling is employed. A description of the constitutive model, the gradient enhancement extension and the implementation into an open-access Abaqus user element subroutine is provided. Parametric studies on the robustness of the approach are demonstrated using two benchmark examples. Finally, the presented approach is used to simulate a detailed femur model.

**Keywords** Bone remodelling · Finite element method · Checkerboarding · Gradient enhancement · Abaqus user element

# 1 Introduction

## 1.1 Background

Bones are metabolically active tissue undergoing constant renewal processes triggered by external and internal factors like proteins, hormonal or mechanical stimuli. In particular, there is an important sensitivity to the mechanical loading conditions which leads to the adaptation of the bone structure, first stated in Wolff's law [1–3]. This continuous adaptation process is referred to as bone remodelling. Bones form optimal structures, which, with minimal material effort, form an optimal load-bearing structure. For instance, active individuals typically exhibit stronger bones, whereas astronauts tend to experience a loss of bone mass density (BMD). At the cellular scale, three different cell types are primarily involved in this process. Osteocytes detect the need for bone repair and renewal. Subsequently, these cells transmit signals to osteoblasts, responsible for removing defective bone tissue, and osteoblasts, responsible for the formation of new bone tissue [2, 4–6].

First attempts to simulate bone remodelling through finite element method (FEM) simulations trace back to the early 90's when predominantly phenomenological models have

Fynn Bensel and Marlis Reiber have contributed equally to this work.

<figure><img src="image_2.png" /></figure>

Fynn Bensel
fynn.bensel@ibnm.uni-hannover.de

Marlis Reiber
marlis.reiber@ibnm.uni-hannover.de

<footnote n="1">¹ Institute of Mechanics and Computational Mechanics (IBNM), Leibniz University Hannover, Appelstraße 9a, Hannover 30167, Germany</footnote>

<footnote n="2">² Institute of Continuum Mechanics (IKM), Leibniz University Hannover, An der Universität 1, Garbsen 30823, Germany</footnote>

<footnote n="3">³ International Research Training Group (IRTG) 2657, Leibniz University Hannover, Appelstraße 11/11a, Hannover 30167, Germany</footnote>

<footnote n="4">⁴ TRR 298: Safety Integrated and Infection Reactive Implants (SIIRI), Hannover Medical School, Carl-Neuberg-Straße 1, Hannover 30625, Germany</footnote>

<footnote n="5">⁵ École normale supérieure Paris-Saclay, Université Paris-Saclay, Avenue des Sciences 4, Gif-sur-Yvette 91190, France</footnote>
been employed, e.g. Carter et al. [7], Beaupré et al. [8], Huiskes et al. [9], Weinans et al. [10], Nackenhorst [11] and many others. These early and purely phenomenological models have already been capable of replicating the typical internal bone structure in response to mechanical loads despite their predominant use of two-dimensional models. To include various effects and dependencies, such as anisotropic behaviour, these models can be adapted further [12, 13]. As computational power increased and more efficient numerical algorithms became available, more sophisticated models have been developed, focusing on smaller length scales, including the renewal of bone tissue down to the cellular level, with the primary objective of gaining a better understanding of the bio-physical processes involved [14–17]. In recent research, multi-scale models have been developed to bridge these scales. Due to the substantial numerical effort of the multi-scale models, their application for patient-specific computations is not feasible yet [18–20].

Nonetheless, already at the phenomenological continuum scale, deficiencies in the stability of the numerical schemes have been observed in the early stages, which lead to the so-called checkerboard phenomenon [21–24]. The phenomenon is evident in some spatial domains where the computed BMD alternates between maximum and minimum values in neighbouring elements instead of exhibiting a smooth representation.

## 1.2 Related research

Different approaches to suppress the formation of checkerboarding in FEM simulations on bone remodelling have been suggested. The approaches range from purely numerical approaches to rigorous mathematical analysis.

### 1.2.1 Filtering

A straightforward approach is to filter the resulting local BMD field by extrapolation and averaging onto the nodes for post-processing [10]. This method corresponds to an image processing filter. Filtering generally does not prevent checkerboarding since the underlying mathematical problem remains ill-posed.

### 1.2.2 Smoothing

Another approach is to apply smoothing in each iteration, e.g. on the mechanical stimulus [21], the strain energy density [24], stress and BMD distribution [25] or via the introduction of spatial influence functions [26]. These so-called node-based smoothing techniques may have the disadvantage of being mesh-dependent or sensitive to integration schemes. Consequently, convergence may not be guaranteed [27]. Calvo-Gallego et al. [27] have recently proposed a novel

non-local smoothing technique to overcome the aforementioned problems.

In general, smoothing methods can be viewed as a remedy for the occurring discontinuities but not as a treatment of their origin. Smoothing does not directly impact the computational cost of the FEM solution, but each iteration step needs to be post-processed, which increases the computational time.

### 1.2.3 Higher-order ansatz functions

Furthermore, higher-order finite elements, e.g. quadratic ansatz functions can be used [28, 29]. This approach has the disadvantage of an increased numerical effort.

### 1.2.4 Nodal BMD formulation

Most models define the BMD as an internal variable so that the balance of mass can be treated locally. In contrast, some approaches suggest to introduce the BMD as a nodal quantity to ensure continuous interpolation. These nodal-based formulations require the simultaneous solution of the balance of mass and momentum, often using a monolithic approach. Further, the different time scales of the change of the displacement and the BMD field need to be addressed. These approaches increase the computational effort and might cause spurious oscillations [21, 28, 30]. Further, Garzón-Alvarado and Linero [31] show that the solution can depend on the mesh and time step size.

### 1.2.5 Mathematical adaptions

The mathematical approaches aim to provide uniqueness and sufficient conditions or supply mathematical extensions to guarantee those [32–34]. Harrigan and Hamilton [35] show that by introducing an additional factor in the evolution equation, the stability and uniqueness of the problem can be restored. This approach needs at least quadratic elements and therefore increases the numerical effort [29].

### 1.2.6 Topology optimisation

The similarity of the modelling of growth to topology optimisation has been acknowledged for a long time [36]. Topology optimisation aims to find the optimal material distribution in a constrained region. Various strategies have been suggested, like the Solid Isotropic Material with Penalisation (SIMP) approach [37]. A new variational growth approach to topology optimisation using regularisation has been introduced by Junker and Hackl [38]. The authors propose a gradient-enhanced strain energy formulation to prevent checkerboarding. Similar techniques have been established in damage mechanics to prevent mesh dependency [39–43].
## 1.3 Scope of the current work

In this work, the method of gradient enhancement is introduced to prevent checkerboarding in bone remodelling simulations. For the FEM simulation the phenomenological bone remodelling theory from Lutz and Nackenhorst [24] is adopted. Results of the latest research justify the application of phenomenological models based on the strain energy density stimulus [3, 44].

The remainder of the article is structured as follows: in Sect. 2, the bone remodelling material model and the gradient enhancement approach are explained. In addition, their numerical treatment and implementation in Abaqus are provided. Furthermore, the occurrence of checkerboarding and its origin are discussed in Sect. 3. Next, in Sect. 4 the results of benchmark calculations are presented, including studies on the newly introduced model parameters and different spatial discretisation levels. Further, results from simulations of a human femur for different initial conditions are shown and discussed. Finally, conclusions are drawn in Sect. 5.

# 2 Methods

## 2.1 Modelling framework

First, the phenomenological constitutive model for simulating bone remodelling at the continuum length scale is outlined. Concerning the thermodynamic consistency of the continuum model within the framework of constitutive modelling of materials it is referred to Lutz (2011) [45] and Nackenhorst (2018) [46]. The following basic assumptions are adopted:

* The small strain theory is assumed since bones do not experience strain of more than 0.4% under usual loading, besides that, its failure strain is about 2.0% [46, 47].

* The entire bone remodelling process occurs over several years. Therefore, the remodelling process is modelled as a quasi-static process [13, 28].

* The temperature is assumed to be constant in the human body. Therefore, the process is modelled as isothermal.

Furthermore, it is assumed that the free energy density $\psi$ depends only on the two internal state variables: the elastic strain $\varepsilon$ and the BMD $\rho$. Thus, the strain energy density is defined as

$$ \Psi(\rho, \varepsilon) = \rho \psi(\rho, \varepsilon). \quad (1) $$

From the Clausius-Duhem inequality, it can be concluded that the Cauchy stress $\sigma$ can be calculated by

$$ \sigma = \frac{\partial \psi}{\partial \varepsilon}. \quad (2) $$

The non-linear constitutive relation between the BMD and Young's modulus $E$ is defined according to Lutz and Nackenhorst [24] as

$$ E = E_0 \left( \frac{\rho}{\rho_0} \right)^2, \quad (3) $$

where $E_0 = 6500 \text{ N/mm}^2$ and $\rho_0 = 1 \text{ g/cm}^3$ are the corresponding reference values and the BMD is bounded by physiological limits of the minimal and maximal BMD. These limits represent the range from the lowest stiffness of the bone marrow over the cancellous or spongy bone up to the highest stiffness at the cortical or compact bone:

$$ \rho \in [\rho_{\min}, \rho_{\max}]. \quad (4) $$

The minimum value is chosen as $\rho_{\min} = 0.001 \text{ g/cm}^3$ to avoid division by zero (see e.g. Eq. (7)) and prevent ill-conditioning of the stiffness matrix. Further, the maximum is specified as $\rho_{\max} = 2 \text{ g/cm}^3$, which is a generally accepted value for the maximum apparent density of compact bone [48].

The relation in Eq. (3) shows good agreement with experimental data [46]. This relation is inserted into the generalised Hooke's law

$$ \sigma = C \cdot \varepsilon, \quad (5) $$

with the linear elastic material tensor $C^{\text{LE}}$, the following constitutive relation is obtained

$$ \sigma = \left( \frac{\rho}{\rho_0} \right)^2 C^{\text{LE}} \cdot \varepsilon. \quad (6) $$

Thus, the free energy density is obtained

$$ \begin{aligned} \psi &= \left( \frac{\rho}{\rho_0} \right)^2 \psi^{\text{LE}} \\ &= \frac{1}{\rho} \left( \frac{\rho}{\rho_0} \right)^2 \left[ \frac{\lambda}{2} \text{tr}(\varepsilon)^2 + \mu \text{tr}(\varepsilon^2) \right], \end{aligned} \quad (7) $$

where $\lambda$ and $\mu$ are the Lamé parameters.

From the balance of mass, omitting mass fluxes, the evolution equation for the BMD follows as

$$ \frac{\partial \rho}{\partial t} = \dot{\rho}. \quad (8) $$

Hereby, $t$ could also denote the process time in the scope of a quasi-static simulation. The mass source $\dot{\rho}$ and consequently, the evolution equation for the BMD is defined
according to the strain energy density driven bone remodelling formulation by Beaupré et al. [8] using a first-order approach

$$ \dot{\rho} = c (\Psi - \Psi^{\text{ref}}) = c (\varrho \psi - \Psi^{\text{ref}}). \quad (9) $$

Herein, c is a model parameter that describes the speed of the remodelling process and $\Psi^{\text{ref}}$ is a physiological target value.

## 2.2 Gradient enhancement

The method of a gradient-enhanced strain energy density as used in damage mechanics and topology optimisation is adopted here [38, 39]. Therefore, the strain energy density functional in Eq. (1) is extended by two additional terms

$$ \Psi(\boldsymbol{\epsilon}, \varrho, \phi) = \varrho \psi = \varrho \psi_{\text{mech}}(\boldsymbol{\epsilon}, \varrho) + \frac{\alpha}{2}(\phi - \varrho)^2 + \frac{\beta}{2} |\nabla\phi|^2, \quad (10) $$

where $\alpha$ and $\beta$ are two new model parameters and $\phi$ the nodal representation of the BMD field. This extended strain energy density is used in Eq. (9). Hence, the gradient enhancement directly impacts the evolution equation. A direct penalisation of the local BMD $\varrho$ is avoided since this would impose a $C^1$-continuity requirement on the displacement field $\mathbf{u}$ [39]. Therefore, a global BMD $\phi$ is introduced, which is coupled to the local BMD $\varrho$ via the first additional term. The coupling term forces both BMD fields to coincide, with $\alpha$ controlling this alignment. The second part penalises gradients in $\phi$, with $\beta$ controlling the degree of penalisation.

The global BMD formulation $\phi$ is approximated by $C^0$-continuous ansatz functions. For simplicity, the same ansatz functions are used to discretise the displacement field $\mathbf{u}$ and the BMD field $\phi$.

Equation (10) can be summarised as

$$ \Psi(\boldsymbol{\epsilon}, \varrho, \phi) = \Psi_{\text{mech}}(\boldsymbol{\epsilon}, \varrho) + \Psi_{\text{GE}}(\varrho, \phi), \quad (11) $$

here, $\Psi_{\text{mech}}$ represents the constitutive model and $\Psi_{\text{GE}}$ the gradient enhancement extension.

## 2.3 Numerical treatment

The total energy $\Pi$ is defined as

$$ \Pi := \int_{\Omega} \Psi(\boldsymbol{\epsilon}, \varrho, \phi) \, dV - \int_{\Omega} \mathbf{b} \cdot \mathbf{u} \, dV \\ - \int_{\partial\Omega} \mathbf{t} \cdot \mathbf{u} \, dA \rightarrow \min_{\mathbf{u}, \phi} \quad (12) $$

and has to be minimised with respect to the displacement $\mathbf{u}$ and the global BMD $\phi$. The body forces are represented by $\mathbf{b}$ and $\mathbf{t}$ is the surface traction. This minimisation leads to two variational equations which have to be solved simultaneously

$$ \int_{\Omega} \frac{\partial \Psi}{\partial \boldsymbol{\epsilon}} : \delta \boldsymbol{\epsilon} \, dV - \int_{\Omega} \mathbf{b} \cdot \delta \mathbf{u} \, dV \\ - \int_{\partial \Omega} \mathbf{t} \cdot \delta \mathbf{u} \, dA \stackrel{!}{=} 0 \quad \forall \delta \mathbf{u} \quad (13) $$

and

$$ \int_{\Omega} \alpha (\phi - \varrho) \delta\phi \, dV \\ + \int_{\Omega} \beta \nabla\phi \nabla\delta\phi \, dV \stackrel{!}{=} 0 \quad \forall \delta\phi. \quad (14) $$

The variational equations are solved using FEM. Therefore, the equations are discretised using

$$ \mathbf{u}(\mathbf{x}) = \mathbf{N}_{\mathbf{u}} \hat{\mathbf{u}}, \quad (15) $$

$$ \boldsymbol{\epsilon}(\mathbf{x}) = \mathbf{B} \hat{\mathbf{u}}, \quad (16) $$

$$ \phi(\mathbf{x}) = \mathbf{N}_{\phi} \hat{\phi}, \quad (17) $$

and

$$ \nabla\phi(\mathbf{x}) = \nabla\mathbf{N}_{\phi}\hat{\phi} \quad (18) $$

where $\mathbf{N}_{\mathbf{u}}$ and $\mathbf{N}_{\phi}$ are linear shape functions for the respective solution fields, $\mathbf{B}$ is the shape function derivative matrix and $\hat{\mathbf{u}}$ and $\hat{\phi}$ are the nodal values. The variations of the displacement $\delta\mathbf{u}$, the BMD $\delta\phi$ as well as the strain $\delta\boldsymbol{\epsilon}$ and BMD gradient $\nabla\delta\phi$ can be described by

$$ \delta\mathbf{u} = \mathbf{N}_{\mathbf{u}} \delta\hat{\mathbf{u}}, \quad (19) $$

$$ \delta\phi = \mathbf{N}_{\phi} \delta\hat{\phi}, \quad (20) $$

$$ \delta\boldsymbol{\epsilon} = \mathbf{B} \delta\hat{\mathbf{u}} \quad (21) $$

and

$$ \nabla\delta\phi = \nabla\mathbf{N}_{\phi}\delta\hat{\phi}. \quad (22) $$

Assuming zero body forces $\mathbf{b}$, the residual of Eq. (13) can be expressed as

$$ \mathbf{R}_{\mathbf{u}} = \int_{\Omega} \mathbf{B}^{T} \boldsymbol{\sigma} \, dV - \int_{\partial\Omega} \mathbf{N}_{\mathbf{u}}^{T} \mathbf{t} \, dA = \mathbf{0} \quad (23) $$

and of Eq. (14) as

$$ \mathbf{R}_{\phi} = \int_{\Omega} \alpha(\mathbf{N}_{\phi} \hat{\phi} - \varrho)\mathbf{N}_{\phi} \, dV $$
$$+ \int_{\Omega} \beta \nabla \mathbf{N}_{\phi}^{T} \nabla \mathbf{N}_{\phi} \hat{\boldsymbol{\phi}} \, dV = \mathbf{0}. \qquad (24)$$

This system is solved monolithically. Its incremental form is expressed as

$$\begin{bmatrix} \mathbf{R}_u \\ \mathbf{R}_\phi \end{bmatrix}^{(i+1)} = \mathbf{R}^{(i+1)} \\ = \mathbf{R}^{(i)} + \mathbf{A}^{(i+1)} \begin{bmatrix} \Delta \hat{\mathbf{u}} \\ \Delta \hat{\boldsymbol{\phi}} \end{bmatrix} = \mathbf{0} \qquad (25)$$

where (i) indicates the process time step. The system matrix $\mathbf{A}^{(i+1)}$ is defined as

$$\mathbf{A}^{(i+1)} = \begin{bmatrix} \frac{\partial \mathbf{R}_u^{(i+1)}}{\partial \hat{\mathbf{u}}^{(i+1)}} & \frac{\partial \mathbf{R}_u^{(i+1)}}{\partial \hat{\boldsymbol{\phi}}^{(i+1)}} \\ \frac{\partial \mathbf{R}_\phi^{(i+1)}}{\partial \hat{\mathbf{u}}^{(i+1)}} & \frac{\partial \mathbf{R}_\phi^{(i+1)}}{\partial \hat{\boldsymbol{\phi}}^{(i+1)}} \end{bmatrix}. \qquad (26)$$

The entries from $\mathbf{A}^{(i+1)}$ are calculated as follows

$$\begin{align} \frac{\partial \mathbf{R}_u^{(i+1)}}{\partial \hat{\mathbf{u}}^{(i+1)}} &= \int_{\Omega} \mathbf{B}^T \frac{\partial \mathbb{C}^{(i+1)} \cdot \boldsymbol{\epsilon}^{(i+1)}}{\partial \boldsymbol{\epsilon}^{(i+1)}} \frac{\partial \boldsymbol{\epsilon}^{(i+1)}}{\partial \hat{\mathbf{u}}^{(i+1)}} \, dV \nonumber \\ &= \int_{\Omega} \mathbf{B}^T \left( \mathbb{C}^{(i+1)} + \frac{\partial \mathbb{C}^{(i+1)}}{\partial \mathcal{Q}^{(i+1)}} \right. \nonumber \\ & \qquad \left. \times \frac{\partial \mathcal{Q}^{(i+1)}}{\partial \boldsymbol{\epsilon}^{(i+1)}} \boldsymbol{\epsilon}^{(i+1)} \right) \mathbf{B} \, dV \nonumber \\ &= \int_{\Omega} \mathbf{B}^T \mathbb{C}_{\text{alg}}^{(i+1)} \mathbf{B} \, dV \tag{27} \end{align}$$

and

$$\frac{\partial \mathbf{R}_u^{(i+1)}}{\partial \hat{\boldsymbol{\phi}}^{(i+1)}} = \int_{\Omega} \mathbf{B}^T \frac{\partial \mathbb{C}^{(i+1)}}{\partial \mathcal{Q}^{(i+1)}} \boldsymbol{\epsilon}^{(i+1)} \nonumber \\ \times \frac{\partial \mathcal{Q}^{(i+1)}}{\partial \hat{\boldsymbol{\phi}}^{(i+1)}} \, dV \qquad (28)$$

and

$$\frac{\partial \mathbf{R}_{\phi}^{(i+1)}}{\partial \hat{\mathbf{u}}^{(i+1)}} = \int_{\Omega} -\alpha \mathbf{N}_{\phi} \otimes \frac{\partial \mathcal{Q}^{(i+1)}}{\partial \boldsymbol{\epsilon}^{(i+1)}} \mathbf{B} \, dV \qquad (29)$$

and

$$\frac{\partial \mathbf{R}_{\phi}^{(i+1)}}{\partial \hat{\boldsymbol{\phi}}^{(i+1)}} = \int_{\Omega} \alpha \left( \mathbf{N}_{\phi} - \frac{\partial \mathcal{Q}^{(i+1)}}{\partial \hat{\boldsymbol{\phi}}^{(i+1)}} \right) \otimes \mathbf{N}_{\phi} \, dV \\ + \int_{\Omega} \beta \nabla \mathbf{N}_{\phi}^{T} \nabla \mathbf{N}_{\phi} \, dV. \qquad (30)$$

The internal variable update to compute the local BMD field is performed using an implicit Euler scheme, which

yields

$$\begin{align} \varrho^{(i+1)} &= \varrho^{(i)} + \Delta\varrho \\ &= \varrho^{(i)} + \Delta t \cdot \dot{\rho}(\boldsymbol{\epsilon}^{(i+1)}, \hat{\boldsymbol{\phi}}^{(i+1)}). \tag{31} \end{align}$$

By that, the introduced derivatives are obtained by

$$\frac{\partial \mathbb{C}^{(i+1)}}{\partial \varrho^{(i+1)}} = 2 \frac{\varrho^{(i+1)}}{\varrho_0^2} \mathbb{C}^{\text{LE}}, \qquad (32)$$

and

$$\begin{align} \frac{\partial \varrho^{(i+1)}}{\partial \boldsymbol{\epsilon}^{(i+1)}} &= \frac{\partial \Delta \varrho}{\partial \boldsymbol{\epsilon}^{(i+1)}} \nonumber \\ &= c \Delta t \left( \frac{\varrho^{(i)}}{\varrho_0} \right)^2 \cdot \mathbb{C}^{\text{LE}} \boldsymbol{\epsilon}^{(i+1)} \tag{33} \end{align}$$

and

$$\begin{align} \frac{\partial \varrho^{(i+1)}}{\partial \hat{\boldsymbol{\phi}}^{(i+1)}} &= \frac{\partial \Delta \varrho^{(i+1)}}{\partial \hat{\boldsymbol{\phi}}^{(i+1)}} \nonumber \\ &= c \Delta t \left( \alpha \mathbf{N}_{\phi}^{T} \left( \mathbf{N}_{\phi} \hat{\boldsymbol{\phi}}^{(i+1)} - \varrho^{(i)} \right) \right. \nonumber \\ & \qquad \left. + \beta \nabla \mathbf{N}_{\phi}^{T} \nabla \mathbf{N}_{\phi} \hat{\boldsymbol{\phi}}^{(i+1)} \right). \tag{34} \end{align}$$

### 2.3.1 Implementation in Abaqus

The commercial FEM software Abaqus (Abaqus 2017, Dassault Systèmes, Vélizy-Villacoublay, France) is used for the numerical solution.

A user element (UEL) subroutine is used to implement the gradient-enhanced bone remodelling. The basis is the Abaqus element C3D4, a four-node linear tetrahedral element. Tetrahedral elements are chosen because of their ability to mesh complicated geometries automatically. Nevertheless, a vast number of elements are needed to represent the geometry of bones, for example, the femur. The fine mesh is needed for a high resolution of the geometry and the BMD field rather than the approximation of the displacement gradients. As the displacement gradients are well-approximated by the fine mesh, the use of quadratic tetrahedral elements is not necessary.

For the gradient enhancement, an additional degree of freedom for the global BMD $\phi$ is added to the three displacement degrees of freedom. The local BMD field $\varrho$ is described as an internal state-dependent variable (SDV). The UEL is depicted as a pseudocode in Appendix A.

Unfortunately, Abaqus does not support the direct visualisation of user elements. This problem can be circumvented by creating a copy of each element. These additional dummy elements are superimposed on the user elements. By using
so-called “common blocks”, the computed results can be transferred internally from the UEL elements to the dummy elements and visualised using a user material subroutine (UMAT) [49, 50].

Alternatively, a thermo-mechanical coupled analysis could be used in combination with a User Material Heat Transfer (UMATHT) for the gradient enhancement [51]. Nonetheless, the UEL approach provides greater flexibility to the user.

## 2.4 Virtual X-ray imaging

X-Ray images are a common tool for the visual assessment of bone quality. To bridge the gap between numerical investigations and medical practice, the computational results are post-processed via virtual X-Ray imaging.

Numerically, a ray r can be fully characterised by its origin o, direction d and a strictly positive factor λ

$$r(\lambda) = o + d \cdot \lambda, \quad \lambda \ge 0. \qquad (35)$$

X-Ray images measure the difference in the ray intensity between the emission $I_0(E_r)$ and detection plate $I(E_r)$. The intensity of the ray can be described by the Beer-Lambert law as

$$I(E_r) = I_0(E_r) e^{-\int_{\mathbf{r}} \alpha_r(E_r, \varrho) \, d\mathbf{r}} \qquad (36)$$

where $\alpha_r$ is the linear attenuation coefficient of the crossed material and $E_r$ is the energy of the ray.

Introducing the mass attenuation coefficient $\mu_r$ as

$$\mu_r(E_r) = \frac{\alpha_r(E_r, \varrho(\mathbf{x}))}{\varrho(\mathbf{x})}, \qquad (37)$$

the energy decay can be computed by

$$\frac{I(E_r)}{I_0(E_r)} = e^{-\mu_r(E_r) \int_{\mathbf{r}} \varrho(\mathbf{x}) d\mathbf{r}}. \qquad (38)$$

This formulation reduces the task to integrating the BMD alongside the rays.

The best practice is to define a box enclosing the full femur model. Afterwards, the opposite sides are defined as an emission and a detection plate, which also define the size of the final image. One ray per pixel needs to be sampled alongside the normal direction of the emission plate [52–54]. The factor can be normalised to $\lambda \in [0, 1]$, such that $r(\lambda = 0)$ is located on the emission and $r(\lambda = 1)$ on the detection plate. This emission is depicted in Fig. 1.

<figure><img src="image_3.png" /><figcaption>Fig. 1 Emission of X-rays onto a FEM model</figcaption></figure>

## 3 Checkerboarding

As mentioned above, bone remodelling simulations are prone to numerical instabilities, the checkerboard phenomenon, which represents a non-physical solution. For example see Figs. 3a or 4a.

### 3.1 Enforcement of physiological bounds

FEM, in its essence, is a minimisation of the total energy. To ensure the uniqueness of the solution the problem must be convex. The convexity of a function or functional can be proven by

$$f(\kappa x + (1 - \kappa)y) \leq \kappa f(x) + (1 - \kappa)f(y) \qquad (39)$$

for $\kappa \in [0, 1]$ and all possible inputs or input functions $x, y \in V$, where $V$ denotes an appropriate input or function space, respectively.

The underlying system of partial differential equations can be summarised as

$$\left\{ \begin{array}{ll} \nabla \cdot \boldsymbol{\sigma}(\boldsymbol{\varepsilon}, \varrho) = 0 & \text{in } \Omega \times \mathbb{T} \\ \text{with } \boldsymbol{\varepsilon} = \nabla^s \mathbf{u} & \\ \text{and } \boldsymbol{\sigma} = \mathbb{C}(\varrho) \cdots \boldsymbol{\varepsilon} & \\ \frac{\partial \varrho}{\partial t} - \dot{\rho} = 0 & \text{in } \Omega \times \mathbb{T} \end{array} \right. . \qquad (40)$$

Here, the first and last equation denote the balance of momentum for the static case under the omission of body forces and the balance of mass without mass fluxes, respectively. Further, $\nabla^s$ denotes the symmetric gradient, $\Omega$ the spatial domain and $\mathbb{T}$ the temporal or process time domain. The system needs to be complemented by appropriate boundary and initial conditions [28, 29, 55]

$$\left\{ \begin{array}{ll} \mathbf{u} = \mathbf{u}_D & \text{in } \partial\Omega_D \times \mathbb{T} \\ \mathbf{n} \cdot \boldsymbol{\sigma} = \mathbf{t} & \text{in } \partial\Omega_N \times \mathbb{T} \\ \varrho^{(0)} = \varrho_{\text{init}} & \text{in } \Omega \times 0 \end{array} \right. . \qquad (41)$$
Here, $\partial\Omega_D$ denotes the Dirichlet boundary, $\partial\Omega_N$ denotes the Neumann boundary, **t** and **n** the surface traction and normal direction, respectively. Further, it must hold that $\partial\Omega_D \cup \partial\Omega_N = \partial\Omega$ and $\partial\Omega_D \cap \partial\Omega_N = \emptyset$.

A difficulty arises from the physiological bounds for the BMD, together with the interdependency of the BMD $\varrho$ and the solution for the displacements **u**, see Eq. (9). Equations (40) and (41) are constrained by

$$ \varrho = \begin{cases} \varrho_{\max} & \text{if } \varrho^{(i)} + \Delta\varrho > \varrho_{\max} \\ \varrho_{\min} & \text{if } \varrho^{(i)} + \Delta\varrho < \varrho_{\min} \\ \varrho^{(i)} + \Delta\varrho & \text{else} \end{cases}. \quad (42) $$

A method to directly incorporate the inequalities in Eq. (42) is to restrict the admissible function space for the solution such that the following holds

$$ \mathcal{V}_h^\varrho = \{ \mathbf{u} \in \mathcal{V}_h | \varrho_{\min} \le \varrho(\mathbf{u}) \le \varrho_{\max} \}. \quad (43) $$

This restriction reduces the original space $\mathcal{V}_h$ to a convex set $\mathcal{V}_h^\varrho \subset \mathcal{V}_h$, which is often discussed in the context of variational inequalities [56]. These restricted function spaces are typically difficult to create and prevent the usage of commercial or open source FEM codes.

The solution must originate from a linear function space to enable the usage of standard finite elements, say $\mathbf{u} \in \mathcal{V}_h$. A typical function space, as in the case of the linear tetrahedral element, would be $\mathcal{V}_h = H^1(\Omega)$. This function space would a priori violate Eq. (42). Therefore, Eq. (31) must be adapted to manually enforce the BMD bounds

$$ \varrho^{(i+1)} = \max \left( \varrho_{\min}, \min \left( \varrho_{\max}, \varrho^{(i)} + \Delta\varrho \right) \right). \quad (44) $$

Such an enforcement can cause numerical instabilities like checkerboarding.

## 3.2 Non-convexity of incremental form

The manual enforcement of the physiological bounds as in Eq. (44), with the dependency of the stiffness on the BMD, can cause the original minimisation problem Eq. (12) to become non-convex. For its convexity, the last two integrals in Eq. (12) can be neglected since they are linear in **u**. Further, $\Pi$ would inherit the convexity from $\Psi$.

The contradiction of Eq. (39) proves the non-convexity of the incremental form. A numerical proof is provided in Appendix B.

## 3.3 Convexity of gradient enhancement

The smoothing effect of the gradient enhancement can be easily shown by the adapted strain energy density functional

<table><caption>**Table 1** Benchmark simulation parameters</caption><thead><tr><th>Parameter</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>E<sub>0</sub></td><td>6500.0</td><td>[N/mm<sup>2</sup>]</td></tr><tr><td>ν</td><td>0.3</td><td>[-]</td></tr><tr><td>ρ<sub>0</sub></td><td>1.0</td><td>[g/cm<sup>3</sup>]</td></tr><tr><td>Ψ<sub>ref</sub></td><td>0.0003</td><td>[N/mm<sup>2</sup>]</td></tr><tr><td>c</td><td>0.02</td><td>[s/m<sup>2</sup>]</td></tr><tr><td>ρ<sub>min</sub></td><td>0.001</td><td>[g/cm<sup>3</sup>]</td></tr><tr><td>ρ<sub>max</sub></td><td>2.0</td><td>[g/cm<sup>3</sup>]</td></tr><tr><td>α</td><td>0.0015</td><td>[m<sup>5</sup>/s<sup>2</sup>kg]</td></tr><tr><td>β</td><td>10<sup>-8</sup></td><td>[m<sup>7</sup>/s<sup>2</sup>kg]</td></tr></tbody></table>

Eq. (10). The expression $|\varrho - \varphi|$ of the coupling term converges against zero. Thus, the coupling term has no influence and can be neglected. But the expression $|\nabla\varphi|$ in the gradient penalisation part is generally non-zero (excluding the trivial solution). Since $|\nabla\varphi|$ is non-zero and monotonically increasing its convexity proof is analogous to $f(x) = x^2$. Hence, the functional in Eq. (10) is introduced, which is convex in its highest gradients and thus regularises the solution of Eq. (12).

# 4 Numerical examples

In order to demonstrate the proposed gradient enhancement for bone remodelling, two benchmark examples have been identified and solved. Afterwards, the approach is applied to a full femur model.

## 4.1 Benchmark problems

### 4.1.1 Plate under compression

The first example is a thin three-dimensional plate (10 cm x 10 cm x 1 cm) with five elements in thickness and 153,164 in total, which is loaded under uniaxial compression (see Fig. 2a). The initial BMD is homogeneous with a value of $\varrho_{init} = 1$ g/cm$^3$. The full set of parameters is summarised in Table 1.

The expected outcome would be a split into two struts with high BMDs which support this system. In the simulations, the load is applied in a first step and held constant in a second step until the BMD converges. Due to the single compression load, the result should have clear and distinct areas of minimum, intermediate and maximum BMD and thus should be very prone to checkerboarding.

In Fig. 3a the result without the gradient enhancement is depicted. The two struts are visible, but checkerboarding occurs in the zones of intermediate densities.
<figure><img src="image_4.png" /><figcaption>(a) Plate under compression</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Plate under shear</figcaption></figure>

<caption>Fig. 2 FEM model for plate under compression and plate under shear</caption>

Applying the gradient enhancement, the result for the BMD in Fig. 3b is obtained. The checkerboarding vanishes and a transition zone with intermediate BMDs can be found between the struts of maximal BMD and the zones of minimal BMD.

### 4.1.2 Plate under shear

The same model of Sect. 4.1.1 is tested again with altered boundary conditions. The plate is clamped on one side and subjected to shear loading as depicted in Fig. 2b. All parameters as well as the total applied load are identical, as summarised in Table 1. Also, for this benchmark problem, the gradient enhancement can prevent the formation of checkerboarding compared to the unregularised solution (see Fig. 4a) and can provide a continuous solution resembling the developing support arch structure with maximum BMD as well as regions of intermediate and minimal BMD (see Fig. 4b).

## 4.2 Influence of characteristic length

To investigate the influence of the new model parameters, $\beta$ is varied. Important to state is that the expression $\sqrt{\frac{\beta}{\alpha}}$ has a length unit. This length is referred to as characteristic length, denoted $h_c$, and characterises the width of the developing transition zones. In Fig. 3b, $h_c$ can be seen in the width of the transition zones towards the top and left edge and in Fig. 4b, it characterises the width of the intermediate BMD zones towards the top-right and bottom-right corner. This description is analogue to phase-field fracture (PFF), where $h_c$ characterises half the width of the fracture, which is the distance from the minimum to the maximum value of the introduced PFF variable (here comparable to the global BMD). The difference arises since the width of a full fracture is zero, but for bone remodelling areas of maximum BMD evolve [56]. Thus, $h_c$ could be used to estimate the gradient enhancement parameters as

$$ \alpha \approx h_c \cdot \Psi^{\text{ref}}, \\ \sqrt{\frac{\beta}{\alpha}} \approx h_c. \qquad (45) $$

The characteristic length scale poses a requirement on the mesh since the elements need to be small enough to resolve areas of intermediate BMD, which are governed by $h_c$. The minimum element length should be smaller than the characteristic length.

In Fig. 3b the result for $h_c \approx 3$ mm is shown. The influence of $h_c$ is depicted in Fig. 5 for $h_c^* = 0.5 h_c$ and $h_c^* = 2 h_c$. This adaption has been achieved by maintaining a constant $\alpha$ while varying $\beta$ between $\frac{1}{4}$ and 4 times its initial value. The developing structure with maximum BMD is similar. Increasing $h_c^*$, the gradients are stronger penalised and the transition zone widens. In Fig. 5b the gradient penalisation is too strong and areas of maximum BMD at some parts are smeared. In Fig. 5a $h_c$ is approximately the characteristic element length.

## 4.3 Spatial discretisation

In Fig. 3b the result for an original element size 0.2 cm is shown. Furthermore, the problem is computed with two additional spatial discretisations, one coarser with 0.25 cm and one finer with 0.1 cm element size. All simulations have been performed with the same parameters; see Table 1. The results of this mesh study are shown in Fig. 6. The depicted results are similar for regions of maximum, minimum BMD and transition zones. For the coarsest mesh, slight differences are observable. The differences vanish with increasing discretisation level.
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Results of bone remodelling simulation for plate under shear loading with and without gradient-enhanced free energy</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 5** Results for final BMD φ distribution in g/cm³ after bone remodelling simulation with different characteristic length h<sub>c</sub></figcaption></figure>

<figure><img src="image_8.png" /></figure>

For the convergence study, a reference mesh with 0.05 cm element size is calculated. In Fig. 7 the convergence of the error measure e of the final BMD solution, introduced as

$$e = \frac{\sqrt{\sum_{j=1}^{N_{\text{nodes}}} (\phi^{(j)} - \phi_{\text{ref}}^{(j)})^2}}{N_{\text{nodes}} \cdot \frac{1}{V} \int_{\Omega} \phi_{\text{ref}} \, dV}, \quad (46)$$

where $N_{\text{nodes}}$ is the number of nodes and V the total femur volume, is shown. The error decreases exponentially with increasing degrees of freedom. The mesh study shows that the solution is independent from the spatial discretisation mesh.

Additionally, the problem is solved using quadratic elements without the gradient enhancement. Therefore, the bone remodelling has been implemented using a UMAT. The results are shown in Fig. 8. The expected support structure is developed and no checkerboarding is apparent. However, the width of the intermediate density areas decreases and in some parts nearly vanishes. The solution shows over- and under-shoots of the physiological bounds, particularly in domains with high gradients of the BMD field. The over- and under-shoots can be attributed to the enforcement in Eq. (44) being applied at the integration point level. As in Fig. 8, averaged element values should be depicted, here an obvious post-processing issue is visible.

The resolution of the mesh, or rather the element size, is governed by the application to bones, since the mesh needs to
<figure><img src="image_9.png" /><figcaption>**Fig. 6** Results of bone remodelling simulation for plate under compression with different mesh discretisation levels (BMD $\phi$ in g/cm$^3$)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7** Convergence of error measure $e$ for different mesh fidelities</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8** Results for plate under compression with quadratic elements without gradient enhancement (BMD in g/cm$^3$). Depicted are constant element values</figcaption></figure>

be fine enough to accurately depict the patient's bone and its BMD distribution. Therefore, the number of elements is the same as for the reference shown in Sect. 4.1.1. The usage of linear elements with gradient enhancement reduces the total number of degrees of freedom by a factor of five. In contrast the computing time is reduced by a factor of four, which is attributed to the UMAT that is required for the visualisation of the UEL results.

## 4.4 Femur

The bone remodelling with the gradient enhancement is applied to a human femur model consisting of 406,905 elements with 304,596 degrees of freedom. The boundary conditions are adopted from Lutz and Nackenhorst [24]. These boundary conditions represent a clamping at the bottom and loading by the joint force, which is the main loading, but also include nine muscle forces (gluteus maximus, medius and minimus, vastus lateralis, intermedialis and medialis, psoas major, biceps femoris, adductor longus), resulting in 275 individually distributed nodal forces, see Fig. 9. The loads have been calculated as statically equivalent loads corresponding to the measured BMD distribution mapped from CT-data to the FEM model [24, 45].

To analyse the uniqueness of the solution, different initial homogeneous BMD distributions, namely 0.2 g/cm$^3$, 1.0 g/cm$^3$, 1.8 g/cm$^3$ and a pseudo random uniform distribution with values between 0.2 g/cm$^3$ and 1.8 g/cm$^3$, are chosen. However, in practice, the starting point is always the mapped BMD field. In a first step, the load is applied and afterwards held constant. During the second step, 1000 quasi-static iterations are performed to analyse the convergence of the BMD distribution. The used material parameters can be found in Table 2.

In Fig. 10, the convergence of the absolute variation in BMD and the variation of the total mass are depicted for the different initial conditions. For all initial conditions, the average variation of the BMD converges to zero (see Fig. 10a). The convergence curves are different for the initial conditions at the beginning, but, they all follow the same convergence rate in the end. The total mass converges to the same value for all initial conditions (see Fig. 10b). The value is in the magnitude of the weight of a human femur [57]. For the medium starting value, the mass decreases and then increases again. First, the outer boundary with maximum BMD forms and then the central elements decrease their BMD to the minimal value, leading to a decrease in total mass. Then, an increase is observed again, corresponding to the decrease of the tran-
<figure><img src="image_12.png" /><figcaption>**Fig. 9** FEM femur model with mesh and loads with close up of the femoral head</figcaption></figure>

<table><caption>**Table 2** Femur simulation parameters</caption><thead><tr><td>Parameter</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>E<sub>0</sub></td><td>6500.0</td><td>[N/mm<sup>2</sup>]</td></tr><tr><td>ν</td><td>0.3</td><td>[-]</td></tr><tr><td>ρ<sub>0</sub></td><td>1.0</td><td>[g/cm<sup>3</sup>]</td></tr><tr><td>Ψ<sub>ref</sub></td><td>0.002</td><td>[N/mm<sup>2</sup>]</td></tr><tr><td>c</td><td>0.01</td><td>[s/m<sup>2</sup>]</td></tr><tr><td>ρ<sub>min</sub></td><td>0.001</td><td>[g/cm<sup>3</sup>]</td></tr><tr><td>ρ<sub>max</sub></td><td>2.0</td><td>[g/cm<sup>3</sup>]</td></tr><tr><td>α</td><td>0.01</td><td>[m<sup>5</sup>/s<sup>2</sup>kg]</td></tr><tr><td>β</td><td>10<sup>-8</sup></td><td>[m<sup>7</sup>/s<sup>2</sup>kg]</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Fig. 10** Convergence plots of the average of the absolute variation of BMD values (a) and the total mass (b) for different initial conditions $\rho_{init}$</figcaption></figure>

sition zone. As expected, the mass decreases for the maximal starting value and vice versa for the minimal starting value.

In Fig. 11, the final results from the bone remodelling process of a femur for different initial conditions are shown. The tubular bone structure is clearly visible for all initial conditions, with maximal BMD on the outside and minimal BMD on the inside. In the femoral head, slight influences can be observed due to the joint load distribution on the mesh nodes. To quantitatively compare the final BMD distributions, the error measure $e$ in Eq. (46) is used. The errors for the final BMD distributions are summarised in Table 3. The differences between the homogeneous distributions are between $0.22 \times 10^{-3}$ and $0.46 \times 10^{-3}$. The largest difference can be found for the random distribution, which is attributed to the local differences in the femoral head.

The shown results illustrate that the mechanical stimulus is of primal influence for forming the internal bone struc-

ture. Furthermore, the results justify the applicability of the phenomenological continuum constitutive theory, which has been previously demonstrated [45]. The proposed method can efficiently and accurately compute a BMD distribution that depicts the bone marrow as well as the cortical and trabecular compartments, starting from a non-physiological BMD distribution. The final distribution appears insensitive to the initial conditions. From Fig. 10, it is evident that a good accuracy is obtained already for 110 iterations. For this simulation, the total wall clock time is 1613s (with around 5571s CPU time) on a workstation with 32GB RAM and four CPUs used (Intel Core i7-9700@3.00GHz).

Considering clinical applications, the bone remodelling surrounding artificial implants is of major interest. In practice, as mentioned before, the starting point is a nearly physiological BMD distribution obtained from image data. In this case, only a few equilibrium iterations are needed after the calculation of the statically equivalent loads [45].
<figure><img src="image_14.png" /><figcaption>**Table 3** Error measure $e$ [%c] between the final solutions of the BMD distribution for different initial conditions in g/cm$^3$</figcaption></figure>

<table><thead><tr><th>Initial BMD</th><th>0.2</th><th>1.8</th><th>Random</th></tr></thead><tbody><tr><td>1</td><td>0.34</td><td>0.22</td><td>0.96</td></tr><tr><td>0.2</td><td>-</td><td>0.46</td><td>0.99</td></tr><tr><td>1.8</td><td>-</td><td>-</td><td>0.94</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>**Fig. 12** Virtual X-ray image for FEM result for BMD $\phi$ after bone remodelling simulation for anteroposterior and lateral view</figcaption></figure>

In Fig. 12, the corresponding X-ray image is depicted with a resolution of 150 × 300 pixels with 80 sample points along the ray. The different zones of cancellous and compact bone in the proximal femur are clearly visible.

A further advantage of the global BMD $\phi$ is its immediate definition at the nodes, which simplifies the computation of the sample points for the virtual X-ray emulation.

## 5 Conclusions

In conclusion, an alternative approach to the existing strategies to prevent checkerboarding in bone remodelling simulations has been presented. It has been demonstrated that the approach of gradient enhancement, previously established in damage mechanics and topology optimisation, can be applied to bone remodelling.

Gradient enhancement enables the usage of standard linear finite elements and the computation of unique solutions with respect to the model parameters $\alpha$ and $\beta$. Moreover, the approach is flexible in the definition of the constitutive equations. Thus, it could be easily extended to, e.g. anisotropic material behaviour.

The method can be easily adopted to almost any FEM solver, like the commercial programme Abaqus. An additional degree of freedom needs to be introduced, increasing the computational effort, but the increase is minor compared to other techniques, like quadratic elements. A benefit of the additional degree of freedom is the availability of the BMD directly at the nodes. This eases further post-processing of the results, like e.g. virtual X-Ray imaging.

The presented estimators for the parameters $\alpha$ and $\beta$ are able to provide accurate results and thus offer a ready-to-use application of the associated and open-source UEL subroutine.
# 6 Supplementary information

The commented user element subroutine for Abaqus is published under the Creative Commons Attribution-NonCommercial 3.0 Unported license Ⓒ and available under: https://doi.org/10.25835/th2mqzom.

**Acknowledgements** The funding of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) of the International Research Training Group IRTG 2657 grant 433082294 and the SFB/TRR-298-SIIRI - Project-ID 426335750 is gratefully acknowledged. This work was supported by the compute cluster, which is funded by the Leibniz University of Hannover, the Lower Saxony Ministry of Science and Culture (MWK) and the DFG.

**Author Contributions** The authors FB and MR contributed equally with regard to the writing of the paper and performing the simulations. EF implemented the framework for the X-ray imaging and generated the virtual X-ray images. PJ assisted in the development of the presented method and review of the manuscript. UN provided the idea and assisted with ongoing guidance and review of the paper.

**Funding** Open Access funding enabled and organized by Projekt DEAL. F. Bensel is funded by the DFG through the IRTG 2657 grant 433082294. M. Reiber is funded by the DFG through the SFB/TRR-298-SIIRI - Project-ID 426335750.

**Data availability** An exemplary UEL Abaqus subroutine is available under the aforementioned link.

## Declarations

**Conflict of interest** The authors declare that they have no competing interests of any nature.

**Ethical approval** Not applicable.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## Appendix A Bone remodelling material algorithm

The bone remodelling material algorithm of the UEL is depicted as a pseudocode in Algorithm 1. Before the first increment the SDVs are initialised. The local coordinates $\xi$ and the weights $w$ are extracted. Next, the shape functions $N$ and the local derivatives are calculated $\partial N/\partial x$. Then, a loop over all Gauss points is executed. The Jacobian $J$ and derivative

tives of the shape functions $\partial N/\partial x$ are calculated. Next, the strains $\boldsymbol{\varepsilon}^{(i)}$, the incremental change for the BMD $\Delta\varrho$ and the stresses $\boldsymbol{\sigma}^{(i)}$ are updated. With that, the submatrices and residuals can be evaluated. Finally, the right-hand side (RHS) and the stiffness matrix (AMATRIX) can be assembled. When convergence is reached, the next increment is started, otherwise the next iteration starts. The internal variables SDVs are only updated at the beginning of the next increment.

## Appendix B Non-convexity proof

This example proves the claim from Sect. 3.2 by contradiction of Eq. (39). All introduced constants as well as $\varrho^{(i)}$ are set to one and $\Delta t = 2$ s. Further assuming two solutions $\mathbf{u}_1^{(i+1)}, \mathbf{u}_2^{(i+1)}$ such that $\boldsymbol{\varepsilon}_1^{(i+1)} = 0.25 \mathbb{I}$ and $\boldsymbol{\varepsilon}_2^{(i+1)} = 0.5 \mathbb{I}$ holds. This yields $\bar{\boldsymbol{\varepsilon}}^{(i+1)} = \kappa \boldsymbol{\varepsilon}_1^{(i+1)} + (1-\kappa) \boldsymbol{\varepsilon}_2^{(i+1)} = 0.375 \mathbb{I}$ with $\kappa = 0.5$ and $\mathbb{I}$ being the identity matrix in $\mathbb{R}^3$. Introducing

$$ \psi^{\text{LE}} = \varrho \psi^{\text{LE}} = \left[ \frac{\lambda}{2} \text{tr}(\boldsymbol{\varepsilon})^2 + \mu \text{tr}(\boldsymbol{\varepsilon}^2) \right], $$

it can be computed that

$$ \varrho^{\text{LE}}(\boldsymbol{\varepsilon}_1^{(i+1)}) = 0.75 \ [\text{N/mm}^2], $$

$$ \varrho^{\text{LE}}(\boldsymbol{\varepsilon}_2^{(i+1)}) = 3.0 \ [\text{N/mm}^2], $$

$$ \varrho^{\text{LE}}(\bar{\boldsymbol{\varepsilon}}^{(i+1)}) = 1.6875 \ [\text{N/mm}^2], $$

substituting into Eq. (44) yields

$$ \varrho^{(i+1)}(\boldsymbol{\varepsilon}_1^{(i+1)}) = 0.5 \ [\text{g/cm}^2], $$

$$ \varrho^{(i+1)}(\boldsymbol{\varepsilon}_2^{(i+1)}) = 2.0 \ [\text{g/cm}^2], $$

$$ \varrho^{(i+1)}(\bar{\boldsymbol{\varepsilon}}^{(i+1)}) = 2.0 \ [\text{g/cm}^2], $$

and evaluating Eq. (1) subjected to Eq. (7) gives

$$ \Psi(\boldsymbol{\varepsilon}_1^{(i+1)})|_{\varrho^{(i+1)}} = 0.1875 \ [\text{N/mm}^2], $$

$$ \Psi(\boldsymbol{\varepsilon}_2^{(i+1)})|_{\varrho^{(i+1)}} = 12.0 \ [\text{N/mm}^2], $$

$$ \Psi(\bar{\boldsymbol{\varepsilon}}^{(i+1)})|_{\varrho^{(i+1)}} = 6.75 \ [\text{N/mm}^2]. $$

Substituting this into Eq. (39) yields

$$ 0.5 \cdot 1.875 + 0.5 \cdot 12.0 = 6.09375 < 6.75 \quad \square $$

and thus proves the non-convexity of the incremental form with enforced limits on the BMD.
**Algorithm 1 Pseudocode for bone remodelling UEL**

**1. Initialisation**

Initialise SDVs $\rho_0$
Compute quadrature rule $\xi$, $w$
Get ansatz functions and derivatives $N$, $\partial N/\partial \xi$

**2. Next increment $k += 1$**

Set iteration counter $i = 0$
Update $\rho^{(k)} = \rho^{(k-1)} + \Delta\rho^{(k-1)}$

**3. Next iteration $i += 1$**

**for** loop over Gauss points **do**

Compute **A** and $\partial N/\partial \mathbf{x}$

Update strains $\boldsymbol{\epsilon}$

Calculate $\Delta\rho^{(i)}$

Update stresses $\boldsymbol{\sigma}$

Compute submatrices and residuals $\partial \mathbf{R}_u/\partial \mathbf{u}$, $\partial \mathbf{R}_u/\partial \phi$, $\partial \mathbf{R}_\phi/\partial \mathbf{u}$, $\partial \mathbf{R}_\phi/\partial \phi$, $\mathbf{R}_u$, $\mathbf{R}_\phi$

**end for**

Assemble **AMATRX** and **RHS**

**if** convergence == true **then**

⇒ step 2

**else**

⇒ step 3

**end if**

**References**

1. Wolff J (1892) Das gesetz der transformation der knochen. Hirschwald Verlag, Berlin

2. Robling AG, Castillo AB, Turner CH (2006) Biomechanical and molecular regulation of bone remodeling. Ann Rev Biomed Eng 8:455–498. https://doi.org/10.1146/annurev.bioeng.8.061505.095721

3. Inaba Y, Ike H, Oba M, Saito T (2016) Evaluation of adaptive bone remodeling after total hip arthroplasty using finite element analysis. In: Petrova R (ed) Chap. 9. Perusal of the finite element method. IntechOpen, Rijeka. https://doi.org/10.5772/65031

4. Taber LA (1995) Biomechanics of growth, remodeling, and morphogenesis. Appl Mech Rev 48:487–545. https://doi.org/10.1115/1.3005109

5. Cowin SC (2001) Bone mechanics handbook, 2nd edn. CRC Press, Boca Raton. https://doi.org/10.1201/b14263

6. Bolamperti S, Villa I, Rubinacci A (2022) Bone remodeling: an operational process ensuring survival and bone mechanical competence. Bone Res. https://doi.org/10.1038/s41413-022-00219-8

7. Carter DR, Orr TE, Fyhrie DP (1989) Relationships between loading history and femoral cancellous bone architecture. J Biomech 22(3):231–244. https://doi.org/10.1016/0021-9290(89)90091-2

8. Beaupré GS, Orr TE, Carter DR (1990) An approach for time-dependent bone modeling and remodeling-theoretical development. J Orthop Res 8(5):651–661. https://doi.org/10.1002/jor.1100080506

9. Huiskes R, Weinans H, van Rietbergen B (1992) The relationship between stress shielding and bone resorption around total hip stems and the effects of flexible materials. Clin Orthop Relat Res 274:124–134

10. Weinans H, Huiskes R, Grootenboer HJ (1992) The behavior of adaptive bone-remodeling simulation models. J Biomech 25(12):1425–1441. https://doi.org/10.1016/0021-9290(92)90056-7

11. Nackenhorst U (1997) Numerical simulation of stress stimulated bone remodelling. Tech Mech 17(1):31–40

12. Doblaré M, García JM (2002) Anisotropic bone remodelling model based on a continuum damage-repair theory. J Biomech 35(1):1–17. https://doi.org/10.1016/S0021-9290(01)00178-6

13. Krstin N, Nackenhorst U, Lammering R (2000) Zur konstitutiven Beschreibung des anisotropen beanspruchungsadaptiven Knochenumbaus. Tech Mech 20(1):31–40

14. Lenz C, Nackenhorst U (2004) A numerical approach to mechanosensation of bone tissue based on a micromechanical analysis of a single osteon. In: PAMM: proceedings in applied mathematics and mechanics, vol 4. Wiley, pp 342–343. https://doi.org/10.1002/pamm.200410152

15. Hambli R (2014) Connecting mechanics and bone cell activities in the bone remodeling process: an integrated finite element modeling. Front Bioeng Biotechnol 2:6. https://doi.org/10.3389/fbioe.2014.00006

16. Kardas D, Nackenhorst U, Balzani D (2013) Computational model for the cell-mechanical response of the osteocyte cytoskeleton based on self-stabilizing tensegrity structures. Biomech Model Mechanobiol 12:167–183. https://doi.org/10.1007/s10237-012-0390-y

17. Sapotnick A, Nackenhorst UA (2015) Mechanically stimulated fracture healing model using a finite element framework. In: Lenarz T, Wriggers P (eds) Biomedical technology, vol 74. Springer, Cham, pp 41–53. https://doi.org/10.1007/978-3-319-10981-7_3

18. Webster D, Müller R (2011) In silico models of bone remodeling from macro to nano-from organ to cell. Wiley Interdiscip Rev Syst Biol Med 3(2):241–251. https://doi.org/10.1002/wsbm.115

19. Pastrama M-I, Scheiner S, Pivonka P, Hellmich C (2018) A mathematical multiscale model of bone remodeling, accounting for pore space-specific mechanosensation. Bone 107:208–221. https://doi.org/10.1016/j.bone.2017.11.009

20. García-Aznar JM, Nasello G, Hervas-Raluy S, Ángeles Pérez M, Gómez-Benito MJ (2021) Multiscale modeling of bone tissue mechanobiology. Bone 151:116032. https://doi.org/10.1016/j.bone.2021.116032

21. Jacobs CR, Levenston ME, Beaupré GS, Simo JC, Carter DR (1995) Numerical instabilities in bone remodeling simulations: the advantages of a node-based finite element approach. J Biomech 28(4):449–459. https://doi.org/10.1016/0021-9290(94)00087-K

22. Jantos DR, Hackl K, Junker P (2019) An accurate and fast regularization approach to thermodynamic topology optimization. Int J Numer Methods Eng 117(9):991–1017. https://doi.org/10.1002/nme.5988

23. Dicati GWO, Gubaua JE, Pereira JT (2020) Analysis of the uniqueness and stability of solutions to problems regarding the bone-remodeling process. Med Eng Phys 85:113–122. https://doi.org/10.1016/j.medengphy.2020.10.007

24. Lutz A, Nackenhorst U (2010) Numerical investigations on the biomechanical compatibility of hip-joint endoprostheses. Arch Appl Mech 80(5):503–512. https://doi.org/10.1007/s00419-009-0380-4

25. Gubaua JE, Dicati GWO, da Silva J, do Vale JL, Pereira JT (2022) Techniques for mitigating the checkerboard formation: application in bone remodeling simulations. Med Eng Phys 99:103739. https://doi.org/10.1016/j.medengphy.2021.103739

26. Mullender MG, Huiskes R (1995) Proposal for the regulatory mechanism of Wolff's law. J Orthop Res 13:503

27. Calvo-Gallego JL, Pivonka P, García-Aznar JM, Martínez-Reina J (2021) A novel algorithm to resolve lack of convergence and checkerboard instability in bone adaptation simulations using non-local averaging. Int J Numer Methods Biomed Eng 37(2):3419. https://doi.org/10.1002/cnm.3419

28. Kuhl E, Menzel A, Steinmann P (2003) Computational modeling of growth. Comput Mech 32(1):71–88

29. Bittens M, Nackenhorst U (2023) A fully implicit and thermodynamically consistent finite element framework for bone remodeling
simulations. Comput Mech. https://doi.org/10.1007/s00466-022-02263-x

30. Chen G, Pettet G, Pearcy M, McElwain DLS (2007) Comparison of two numerical approaches for bone remodelling. Med Eng Phys 29(1):134–139. https://doi.org/10.1016/j.medengphy.2005.12.008

31. Garzón-Alvarado DA, Linero D (2012) Comparative analysis of numerical integration schemes of density equation for a computational model of bone remodelling. Comput Methods Biomech Biomed Eng 15(11):1189–1196. https://doi.org/10.1080/10255842.2011.585972

32. Cowin S, Nachlinger RR (1978) Bone remodeling iii: uniqueness and stability in adaptive elasticity theory. J Elast 8:285–295

33. Harrigan TP, Hamilton JJ (1994) Necessary and sufficient conditions for global stability and uniqueness in finite element simulations of adaptive bone remodeling. Int J Solids Struct 31(1):97–107. https://doi.org/10.1016/0020-7683(94)90178-3

34. Harrigan TP, Hamilton JJ (1992) Optimality conditions for finite element simulation of adaptive bone remodeling. Int J Solids Struct 29(23):2897–2906. https://doi.org/10.1016/0020-7683(92)90147-L

35. Harrigan TP, Hamilton JJ (1994) Bone remodeling and structural optimization. J Biomech 27(3):323–328. https://doi.org/10.1016/0021-9290(94)90008-6

36. Bendsøe MP (2003) Aspects of topology optimization and bone-remodelling schemes. Department of Mathematics. Technical University of Denmark. Lyngby, p 16

37. Bendsøe MP (1989) Optimal shape design as a material distribution problem. Struct Optim 1(4):193–202

38. Junker P, Hackl K (2015) A variational growth approach to topology optimization. Struct Multidiscip Optim 52(2):293–304. https://doi.org/10.1007/s00158-015-1241-0

39. Dimitrijevic BJ, Hackl K (2008) A method for gradient enhancement of continuum damage models. Tech Mech Eur J Eng Mech 28(1):43–52

40. Jha NK, Nackenhorst U, Pawar VS, Nadella R, Guruprasad P (2019) On the constitutive modelling of fatigue damage in rubber-like materials. Int J Solids Struct 159:77–89. https://doi.org/10.1016/j.ijsolstr.2018.09.022

41. de Borst R, Mühlhaus H-B (1992) Gradient-dependent plasticity: formulation and algorithmic aspects. Int J Numer Methods Eng 35(3):521–539. https://doi.org/10.1002/nme.1620350307

42. Peerlings RHJ, de Borst R, Brekelmans WAM, de Vree JHP (1996) Gradient enhanced damage for quasi-brittle materials. Int J Numer Methods Eng 39(19):3391–3403. https://doi.org/10.1002/(SICI)1097-0207(19961015)39:19<3391::AID-NME7>3.0.CO; 2-D

43. Polindara C, Waffenschmidt T, Menzel A (2016) Simulation of balloon angioplasty in residually stressed blood vessels-application of a gradient-enhanced fibre damage model. J Biomech 49(12):2341–2348. https://doi.org/10.1016/j.jbiomech.2016.01.037

44. Du J, Li S, Silberschmidt VV (2020) Trabecular bone remodelling: finite-element simulation. Procedia Struct Integrity 28:577–583. https://doi.org/10.1016/j.prostr.2020.10.067. (1st Virtual European Conference on Fracture - VECF1)

45. Lutz A (2011) Ein integrales Modellierungskonzept zur numerischen Simulation der Osseointegration und Langzeitstabilität von Endoprothesen. PhD thesis, Gottfried Wilhelm Leibniz University, Hannover

46. Nackenhorst U (2018) In: Altenbach H, Öchsner A (eds) Modeling of bone adaption processes. Springer, Berlin, pp 1–11. https://doi.org/10.1007/978-3-662-53605-6_33-1

47. Cowin SC, Doty SB (2007) Tissue mechanics. Springer, Heidelberg

48. Martin RB, Burr DB, Sharkey NA, Fyhrie DP (2015) Skeletal tissue mechanics. 2nd edn. Springer, Heidelberg (2015). https://doi.org/10.1007/978-1-4939-3002-9

49. Dassault Systèmes: Abaqus (2016) Online documentation—analysis user's guide. Velizy-Villacoublay, France

50. Jha NK (2016) Visualization of user element in Abaqus. https://doi.org/10.13140/RG.2.2.32755.43043

51. Ostwald R, Kuhl E, Menzel A (2019) On the implementation of finite deformation gradient-enhanced damage models. Comput Mech 64(3):847–877. https://doi.org/10.1007/s00466-019-01684-5

52. Ehlke M, Ramm H, Lamecker H, Hege H-C, Zachow S (2013) Fast generation of virtual x-ray images for reconstruction of 3d anatomy. IEEE Trans Vis Comput Gr 19(12):2673–2682. https://doi.org/10.1109/TVCG.2013.159

53. Li N, Kim S-H, Suh J-H, Cho S-H, Choi J-G, Kim M-H (2007) Virtual X-ray imaging techniques in an immersive casting simulation environment. Nuclear Instrum Methods Phys Res Sect B Beam Interactions Mater Atoms 262(1):143–152. https://doi.org/10.1016/j.nimb.2007.04.262

54. Freud N, Duvauchelle P, Babot D (2003) New developments in virtual x-ray imaging: fast simulation using a deterministic approach. AIP Conf Proc 657(1):553–560. https://doi.org/10.1063/1.1570185

55. Epstein M, Maugin GA (2000) Thermomechanics of volumetric growth in uniform bodies. Int J Plast 16(7):951–978. https://doi.org/10.1016/S0749-6419(99)00081-9

56. Wick T (2020) Multiphysics phase-field fracture. De Gruyter, Berlin. https://doi.org/10.1515/9783110497397

57. Merz AL, Trotter M, Peterson RR et al (1956) Estimation of skeleton weight in the living. Am J Phys Anthropol 14:589–609. https://doi.org/10.1002/ajpa.1330140405

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.