# Analytical beam model revealing bone stresses in femur-implant compound structure

Lukas Pircher, Tilman Grünewald, Helga C. Lichtenegger, Christian Hellmich

► To cite this version:

Lukas Pircher, Tilman Grünewald, Helga C. Lichtenegger, Christian Hellmich. Analytical beam model revealing bone stresses in femur-implant compound structure. Mechanics of Advanced Materials and Structures, 2024, 32 (20), pp.5167-5185. 10.1080/15376494.2024.2419997. hal-05217915

HAL Id: hal-05217915

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

Submitted on 18 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
Analytical beam model revealing bone stresses in femur-implant
compound structure

Lukas Pircher<sup>a</sup>, Tilman A. Grünewald<sup>b</sup>, Helga Lichtenegger<sup>c</sup>, and Christian Hellmich<sup>a</sup>

<sup>a</sup> Institute for Mechanics of Materials and Structures, TU Wien (Vienna University of Technology), Karlsplatz 13/202, 1040 Vienna, Austria, <sup>b</sup> Aix-Marseille Univ, CNRS, Centrale Med, Institut Fresnel, Marseille, France, <sup>c</sup> Institute of Physics and Materials Science, BOKU University, Vienna, Austria

ARTICLE HISTORY

Compiled October 17, 2024

ABSTRACT

We here introduce an analytical model for a preclinical femur-implant compound structure. This model extends standard Euler-Bernoulli-Saint-Venant theory by closed-form expressions for such shear stresses in thin-walled components, which arise from elastic material properties changing at the interfaces between bone and implant materials. Particularly noteworthy are stress singularities occurring at the generators of the cylindrical implants, which intersect the long axis of the bone shaft. This is fully consistent with observed major reorientations of bone fibrils following the implantation event, so that our novel mechanical model clearly indicates a case of micro-mechanobiology.

KEYWORDS

analytical mechanics, composite beam theory, shear stress, local 3D equilibrium,
magnesium implants, bone remodeling

**1. Introduction**

Pre-clinical studies are an indispensable element of biomedical research, and animal
studies remain an indispensable part of preclinical research. However, as ethical con-
cerns call for reduction of animal studies, the latter need to be very carefully designed
[1], together with a maximized degree of evaluation and understanding of the outcome
of corresponding investigations. In this context, mathematical modeling plays an ever
increasing role [2], and the present contribution tackles the clearly needed improvement
of corresponding mathematical models. More specifically, we focus on structural me-
chanics models for bone-implant structures occurring in standard animal models, more
precisely on murine femoral shafts penetrated transcortically by magnesium implants
[3,4], see Figure 1. In this context, mathematical models have typically dealt with the
degradation and corrosion process of the magnesium implants [6-8] — whereby Stefan
problems, when cast into a rigorous mass conservation frame, play a central role [9,10].
However, we here set a different focus, one on the mechanical behavior of the magne-
sium implant and the surrounding bone. This is motivated by a more recent discovery
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Preclinical trial of magnesium implants in murine femora: (a) skeleton of rat with right hindlimp; (b) position of magnesium pin in right femoral shaft, (c) microCT slice depicting implant one month after implantation, image reproduced from [5], under license CC BY 4.0</figcaption></figure>

made possible through synchroton small-angle scattering tomography [5]: significant
changes in the bone microstructure evolve in the weeks and months after implantation;
and this is so presumably because of a dramatically changed mechanical environment.
In more detail, the implant causes the originally dominant orientation direction of
fibrils along the long bone axis to shift towards a vortex-type (or “wrapping”) pattern
around the circumference of the implant [11]. Upon progressive implant resorption,
however, the fibril orientation pattern transforms back towards the original configu-
ration [11]. While noting that the 3D solid Finite Element Method (FEM) [12] is the
current gold standard for mechanical organ-implant simulations [13], the current paper
deals with a different approach, due to two main challenges faced by state-of-the-art
Finite Element models constructed from microComputed Tomograms: (i) they do not
allow for direct study of different design parameters such as implant shape or size, as
any corresponding changes necessitate repeated, generally quite expensive numerical
studies; and (ii) when tested on numerical convergence, large Finite element models
of bony objects have been shown to actually fail [14], so that quantitatively relevant
stress estimations obtained from the FEM cannot be taken for granted, or rather,
remain out of reach in many circumstances. In order to effectively tackle these open
challenges, we here introduce an analytical bone-implant model, along the following
lines: After geometrical representation of the femur-implant structure as two inter-
penetrating cylinders, one hollow and one filled, we recall the standard expressions
for the normal stresses encountered in composite beams with cross-sections reflect-
ing the aforementioned bone-implant compound structure. We continue Section 2 by
re-calling the theory of shear stresses in thin-walled beams undergoing torsion and
bending, before proposing corresponding extensions for composite beams of the afore-
mentioned type, along with new expressions for shear stresses in such beams loaded
by axial forces and bending moments. In Section 3, these expressions are then eval-
uated for practically relevant geometrical and material properties, and illustrated in
terms of corresponding stress distributions. Finally, Section 4 contains a discussion in
view of tomograms elucidating microstructural changes in bone tissue in the vicinity
of the implant, followed by the limitations of the present approach and corresponding
$$ \begin{gathered} R_{bone} = 1.95 \text{ mm } [5] \quad | \quad t_{bone} = 1.00 \text{ mm } [5] \quad | \quad R_{imp} = 0.80 \text{ mm } [5] \\ E_{bone} = 25 \text{ GPa } [14] \quad | \quad E_{imp} = 37 \text{ GPa } [15] \end{gathered} $$

<figure><img src="image_3.png" /><figcaption>Figure 2. Definition of system geometry, by means of bone radius $R_{bone}$, bone thickness $t_{bone}$, and implant radius $R_{imp}$; implant thickness $t_{imp}(x)$ and interface-related polar angle $\gamma_{int}(x)$ depend on axial distance $x$ from the implant axis; positions on the bone profile curve are given in terms of the polar angle $\gamma$ and the arc length $s$</figcaption></figure>

potential future activities.

## 2. Mathematical formulation

### 2.1. Geometrical system, elastic properties, as well as mechanical forces and stresses

We consider two penetrating cylinders, one hollow and one filled, with the first representing a cortical bone shaft and the second representing an implant, see Figures 2 and 3. The axes of these cylinders intersect, while lying perpendicular to each other. This geometrical system is defined through the outer radius of the hollow cylinder, $R_{bone}$, by the thickness of this cylindrical shell, $t_{bone}$, and by the radius of the filled implant cylinder, $R_{imp}$, see Figure 2 and Table 1.

The axis of the first cylinder is regarded as beam axis with orientation $e_x$, and the corresponding composite beam, which consists of two materials, cortical bone and magnesium with elastic properties according to Table 1, is subjected to second-order tensorial Cauchy stress states $\sigma$ of the form [16]

$$ \sigma = \sigma_{xx}e_x \otimes e_x + \sigma_{xy}(e_x \otimes e_y + e_y \otimes e_x) + \sigma_{xz}(e_x \otimes e_z + e_z \otimes e_x) \quad (1) $$
with $\sigma_{xx}$, $\sigma_{xy}$, and $\sigma_{xz}$ as the non-vanishing components of the stress tensor; and with $\mathbf{e}_x$, $\mathbf{e}_y$, and $\mathbf{e}_z$ as the unit vectors forming the orthonormal basis depicted in Figure 3. These stress states are associated to internal forces, namely axial forces $N$, shear forces $S_y$ and $S_z$, torsional moments $M_x$, as well as bending moments $M_y$ and $M_z$. Such internal forces can be obtained from ground reaction forces [17] in combination with an inverse dynamics model for the overall musculoskeletal system [18]. The latter reference provides, at the center of the femoral shaft of a rat, the temporal evolution of the aforementioned internal forces over a single trotting stride. Temporal averaging over this stride yields the following values,

$$ 
\begin{array}{l}
N = -8.03 \, \text{N} \\
M_x = -1.15 \, \text{Nmm} \\
S_y = -1.15 \, \text{N} \\
M_y(x=0) = 12.78 \, \text{Nmm} \\
S_z = -0.16 \, \text{N} \\
M_z(x=0) = -1.20 \, \text{Nmm}
\end{array} 
\quad (2) $$

While the forces and the torsional moments will be considered as constant along the beam axis in the vicinity of implant, the standard equilibrium conditions for slender beams undergoing bending, see e.g. [16,19] and Appendix A,

$$ \frac{dM_y}{dx} = S_z \qquad \frac{dM_z}{dx} = -S_y \quad (3) $$

entail correspondingly linear trends of the bending moments,

$$ M_y(x) = M_y(x=0) + S_z x \qquad M_z(x) = M_y(x=0) - S_y x \quad (4) $$

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Resultant forces acting on bone-implant compound structure (femoral shaft portion in blue, with axis direction $\mathbf{e}_x$; magnesium implant in red, with axis direction $\mathbf{e}_y$); quantified in terms of beam theory-related internal force quantities $N$ (axial force), $M_y$ and $M_z$ (bending moments), $S_y$ and $S_z$ (shear forces), and $M_x$ (torsional moment); $\mathbf{e}_x$ is the direction of the bone axis and $\mathbf{e}_y$ is the direction of the implant axis</figcaption></figure>
## 2.2. Normal stresses according to composite beam theory

Following the standard composite beam theory [20–22], where the axial normal strains $\varepsilon_{xx}$ are linearly distributed throughout the cross sections of the composite beam, the normal stress distributions over the cross sections of the investigated bone-implant structure are linked to the axial force $N$ and the bending moments $M_y$ and $M_z$ through the relation

$$ \sigma_{xx}^{N} = \frac{E(x, s)}{E_{bone}} \frac{N}{A_{bone}(x) + \frac{E_{imp}}{E_{bone}} A_{imp}(x)} \quad (5) $$

$$ \sigma_{xx}^{M_y} = \frac{E(x, s)}{E_{bone}} \frac{M_y(x)z(s)}{I_{y,bone}(x) + \frac{E_{imp}}{E_{bone}} I_{y,imp}(x)} \quad (6) $$

$$ \sigma_{xx}^{M_z} = - \frac{E(x, s)}{E_{bone}} \frac{M_z(x)y(s)}{I_{z,bone}(x) + \frac{E_{imp}}{E_{bone}} I_{z,imp}(x)} \quad (7) $$

In Eqs.(5) to (7), $s$ denotes the arc length labeling positions on the midsurface of the bone shell and its extensions into the implant. It reads mathematically as

$$ s(\gamma) = \gamma \left( R_{bone} - \frac{t_{bone}}{2} \right) \leftrightarrow \gamma(s) = \frac{s}{\left( R_{bone} - \frac{t_{bone}}{2} \right)} \quad (8) $$

with the polar angle $\gamma$ as depicted in Figure 2. The arc length is related to the Cartesian coordinates $y$ and $z$ through

$$ y(s) = \cos[\gamma(s)] \times \left( R_{bone} - \frac{t_{bone}}{2} \right) = \cos\left[\frac{s}{\left(R_{bone} - \frac{t_{bone}}{2}\right)}\right] \times \left( R_{bone} - \frac{t_{bone}}{2} \right) \quad (9) $$

and

$$ z(s) = \sin[\gamma(s)] \times \left( R_{bone} - \frac{t_{bone}}{2} \right) = \sin\left[\frac{s}{\left(R_{bone} - \frac{t_{bone}}{2}\right)}\right] \times \left( R_{bone} - \frac{t_{bone}}{2} \right) \quad (10) $$

Furthermore, in Eqs. (5) to (7), $A_{bone}$ denotes the cross-sectional area of bone; $A_{imp}$ denotes the area of implant replacing former bone; $I_{y,bone}$ and $I_{z,bone}$, respectively, are the second-order area moments around the $y$-axis and the $z$-axis of the bone portion within the cross-section through the bone-implant compound; and $I_{y,imp}$ and $I_{z,imp}$, respectively, are the second-order area moments around the $y$-axis and the $z$-axis of the implant portion within cross-section through the bone-implant compound. We explicitly note that the geometrical properties of the cross-sections, such as areas and area moments, explicitly depend on the coordinate of the beam axis, $x$. In this context, the cross-sectional areas and area moments are quantified in terms of the cross-section-specific polar angle of the bone-implant interface, which we denote as $\gamma_{int}$. It can be given in terms of the section-specific implant thickness $t_{imp}$ and the outer radius of the femoral shaft, $R_{bone}$, through

$$ \gamma_{int}(x) = \arcsin \left( \frac{\frac{t_{imp}(x)}{2}}{R_{bone} - \frac{t_{bone}}{2}} \right) \quad (11) $$
The thickness $t_{imp}$ of the cross-sectional portion made up by the implant can, in turn, be expressed through the axial distance from the cylindrical axes intersection point, denoted as $x$, and the implant radius, via

$$t_{imp}(x) = 2\sqrt{(R_{imp})^2 - x^2} \qquad (12)$$

see Figure 2. The cross sectional properties of the hollow bone cylinder without implant, occurring at $x \le R_{imp}$ and $x \ge R_{imp}$, read as

$$A_{bone} = 2\pi \left[ R_{bone} - \frac{t_{bone}}{2} \right] t_{bone} \qquad (13)$$

$$I_{y,bone} = \frac{1}{4}\pi \left[ (R_{bone})^4 - (R_{bone} - t_{bone})^4 \right] \qquad (14)$$

$$I_{z,bone} = \frac{1}{4}\pi \left[ (R_{bone})^4 - (R_{bone} - t_{bone})^4 \right] \qquad (15)$$

$$A_{imp} = I_{y,imp} = I_{z,imp} = 0 \qquad (16)$$

The cross sectional properties of cross sections comprising both bone and implant, occurring at $-R_{imp} < x < R_{imp}$, read as

$$A_{bone}(x) = (2\pi - 4\gamma_{int}(x)) \left[ R_{bone} - \frac{t_{bone}}{2} \right] t_{bone} \qquad (17)$$

$$I_{y,bone}(x) = \frac{1}{4} \left( \pi - 2\gamma_{int}(x) + \sin[2\gamma_{int}(x)] \right) \times \left[ (R_{bone})^4 - (R_{bone} - t_{bone})^4 \right] \qquad (18)$$

$$I_{z,bone}(x) = \frac{1}{4} \left( \pi - 2\gamma_{int}(x) - 2\cos[\gamma_{int}(x)]\sin[\gamma_{int}(x)] \right) \times \left[ (R_{bone})^4 - (R_{bone} - t_{bone})^4 \right] \qquad (19)$$

$$A_{imp}(x) = 4\gamma_{int}(x) \left[ R_{bone} - \frac{t_{bone}}{2} \right] t_{bone} \qquad (20)$$

$$I_{y,imp}(x) = \frac{1}{2} \left( \gamma_{int}(x) - \cos[\gamma_{int}(x)] \sin[\gamma_{int}(x)] \right) \times \left[ (R_{bone})^4 - (R_{bone} - t_{bone})^4 \right] \qquad (21)$$

$$I_{z,imp}(x) = \frac{1}{2} \left( \gamma_{int}(x) + \cos[\gamma_{int}(x)] \sin[\gamma_{int}(x)] \right) \times \left[ (R_{bone})^4 - (R_{bone} - t_{bone})^4 \right] \qquad (22)$$

The cross-sectional compartments hosting bone- and implant-related axial Young's modulus values, $E_{bone}$ and $E_{imp}$, are given by means of the following function of axial
position x and arc length s

$$E(x, s) = E_{imp} - (E_{imp} - E_{bone}) \times \\ \left\{ H \left[ s - \gamma_{int}(x) \times \left( R_{bone} - \frac{t_{bone}}{2} \right) \right] \right. \\ -H \left[ s - \left[ \pi - \gamma_{int}(x) \right] \times \left( R_{bone} - \frac{t_{bone}}{2} \right) \right] \\ +H \left[ s - \left[ \pi + \gamma_{int}(x) \right] \times \left( R_{bone} - \frac{t_{bone}}{2} \right) \right] \\ \left. -H \left[ s - \left[ 2\pi - \gamma_{int}(x) \right] \times \left( R_{bone} - \frac{t_{bone}}{2} \right) \right] \right\} \quad (23)$$

with H denoting the Heaviside step function; $H(\xi < 0) = 0$ and $H(\xi > 0) = 1$. This function allows for an elegant representation of the jump encountered in the elastic modulus when moving from the bone to the implant portion of the considered compound structure.

## 2.3. Axial gradients of normal stresses occurring in local 3D equilibrium conditions

### 2.3.1. Motivation

As a starting point, we follow the classical beam theory approaches [16,19] where gradients in axial forces arise from shear tractions on the beam surfaces (for the present application, the latter may be realized in terms of point loads associated with muscle attachments), and where shear stresses are due to torsional moments as well as to gradients in bending moments. In the present application, the thickness $t_{bone}$ of the cylindrical shell is always constant, and the axial-circumferential shear stresses $\sigma_{xs}$ are considered as not significantly varying over this thickness, so that $\sigma_{xs}$ is uniquely related to the shear flow $F_{xs}$, via

$$F_{xs} = t_{bone}\sigma_{xs}, \quad (24)$$

whereby the suffix “s” relates to the arc length running along the midsurface of the cylindrical shell of thickness $t_{bone}$, see Figure 2. Such shear stresses may arise from torsional loading, as described further below in Section 2.5, and from gradients in axial normal stresses, which in the present case stem from various sources, as described in Sections 2.3.2 to 2.3.4. The coupling between shear stresses and axial normal stresses results from the following specialized form of the 3D local equilibrium condition [16]

$$\frac{\partial \sigma_{xx}}{\partial x} + \frac{\partial \sigma_{xs}}{\partial s} = 0 \quad (25)$$
where the shear stress component along the cylindrical circle is related to the Cartesian shear stress components via

$$ \sigma_{xy} = -\sigma_{xs} \sin \gamma \qquad (26) $$

$$ \sigma_{xz} = \sigma_{xs} \cos \gamma \qquad (27) $$

and where the gradients of the axial normal stresses will be, in the following, determined from Eqs. (5) to (7). These gradients arise from changes of internal forces, from changes of cross-sectional geometrical properties (areas and area moments), as well as from changes of mechanical properties (values for Young's modulus).

### 2.3.2. Axial gradients of normal stresses due to bending around y-axis

As regards the gradients of axial normal stress due to bending around the y-axis, Eq. (6), together with Eq. (3), implies that

$$ \frac{\partial \sigma_{xx}^{M_y}}{\partial x} = \left( \frac{\partial \sigma_{xx}^{M_y}}{\partial x} \right)^{dM_y} + \left( \frac{\partial \sigma_{xx}^{M_y}}{\partial x} \right)^{dA} + \left( \frac{\partial \sigma_{xx}^{M_y}}{\partial x} \right)^{dE} \qquad (28) $$

where new expressions for axial normal stress gradients arising from a gradient of the bending moment around the y-axis, from gradients of bone- and implant-specific cross-sectional areas and area moments, and from a gradient in the elastic modulus have been introduced. In more detail, in Eq. (28), the normal stress gradients arising from a change in the bending moment are given by the following mathematical expression

$$ \left( \frac{\partial \sigma_{xx}^{M_y}}{\partial x} \right)^{dM_y} = \frac{E(x, s) z(s) S_z}{E_{bone} I_{y,bone}(x) + I_{y,imp}(x) (E_{imp}/E_{bone})} \qquad (29) $$

with $E(x, s)$ according to Eq. (23). Moreover, in Eq. (28), the normal stress gradients arising from changes in the cross-sectional geometrical properties are given by the following mathematical expression

$$ \left( \frac{\partial \sigma_{xx}^{M_y}}{\partial x} \right)^{dA} = - \frac{E(x, s) z(s) M_y(x)}{E_{bone} [I_{y,bone}(x) + I_{y,imp}(x) (E_{imp}/E_{bone})]^2} \times \frac{dI_{y,bone}(x)}{dx} \left[ 1 - \frac{E_{imp}}{E_{bone}} \right] \qquad (30) $$
with the gradient of the second-order area moment reading as

$$ \frac{dI_{y,bone}(x)}{dx} = -\frac{dI_{y,imp}(x)}{dx} \\ = \frac{x}{2} \left[ (R_{bone})^4 - (R_{bone} - t_{bone})^4 \right] \times \\ \times \left( \frac{1 - \cos \left[ 2 \arcsin \left[ \frac{\sqrt{(R_{imp})^2 - x^2}}{R_{bone} - \frac{t_{bone}}{2}} \right] \right]}{R_{bone} - \frac{t_{bone}}{2} \sqrt{(R_{imp})^2 - x^2} \sqrt{1 - \frac{(R_{imp})^2 - x^2}{R_{bone} - \frac{t_{bone}}{2}}}} \right) \quad (31) $$

Finally, in Eq. (28), the normal stress gradients arising from a change in the elastic modulus are given by the following mathematical expression

$$ \left( \frac{\partial \sigma_{xx}^{M_y}}{\partial x} \right)^{dE} = \frac{\partial E(x, s)}{\partial x} \frac{1}{E_{bone}} \frac{z(s) M_y(x)}{I_{y,bone}(x) + \frac{E_{imp}}{E_{bone}} I_{y,imp}(x)} \quad (32) $$

whereby the gradient of the mechanical properties follows from Eq. (23), and consequently, it reads as

$$ \frac{\partial E(x, s)}{\partial x} = \frac{d\gamma_{int}(x)}{dx} \left( R_{bone} - \frac{t_{bone}}{2} \right) (E_{imp} - E_{bone}) \times \\ \left\{ \begin{aligned} & \delta \left[ s - \gamma_{int}(x) \times \left( R_{bone} - \frac{t_{bone}}{2} \right) \right] \\ & + \delta \left[ s - \left[ \pi - \gamma_{int}(x) \right] \times \left( R_{bone} - \frac{t_{bone}}{2} \right) \right] \\ & + \delta \left[ s - \left[ \pi + \gamma_{int}(x) \right] \times \left( R_{bone} - \frac{t_{bone}}{2} \right) \right] \\ & + \delta \left[ s - \left[ 2\pi - \gamma_{int}(x) \right] \times \left( R_{bone} - \frac{t_{bone}}{2} \right) \right] \end{aligned} \right\} \quad (33) $$

In Eq. (33), $\delta$ denotes the Dirac delta function, which is the derivative of the Heaviside function; $\delta(\xi = 0) = +\infty$, $\delta(\xi \neq 0) = 0$ and $\int_{-\infty}^{\infty} \delta(\xi) ds = 1$. Moreover, in Eq.(33), it was also considered that

$$ \frac{d\gamma_{int}(x)}{dx} = - \frac{x}{\sqrt{(R_{imp})^2 - x^2} \sqrt{1 - \frac{(R_{imp})^2 - x^2}{R_{bone} - \frac{t_{bone}}{2}}}} \quad (34) $$

### 2.3.3. Axial gradients of normal stresses due to bending around z-axis

As regards the gradients of axial normal stress due to bending around the z-axis, Eq. (7), together with Eq. (3), implies that

$$ \frac{\partial \sigma_{xx}^{M_z}}{\partial x} = \left( \frac{\partial \sigma_{xx}^{M_z}}{\partial x} \right)^{dM_y} + \left( \frac{\partial \sigma_{xx}^{M_z}}{\partial x} \right)^{dA} + \left( \frac{\partial \sigma_{xx}^{M_z}}{\partial x} \right)^{dE} \quad (35) $$
with new expressions analogous to those introduced below Eq. (28). In more detail, in Eq. (35), the normal stress gradients arising from a change of the bending moment are given by the following mathematical expression

$$ \left( \frac{\partial \sigma_{xx}^{M_z}}{\partial x} \right)^{dM_z} = \frac{E(x, s)}{E_{bone}} \frac{y(s) S_y}{I_{z,bone}(x) + I_{z,imp}(x)(E_{imp}/E_{bone})} \quad (36) $$

with $E(x, s)$ according to Eq. (23). Moreover, in Eq. (35), the normal stress gradients arising from changes in the cross-sectional geometrical properties are given by the following mathematical expression

$$ \left( \frac{\partial \sigma_{xx}^{M_z}}{\partial x} \right)^{dA} = \frac{E(x, s)}{E_{bone}} \frac{y(s) M_z(x)}{[I_{z,bone}(x) + I_{z,imp}(x)(E_{imp}/E_{bone})]^2} \times \frac{dI_{z,bone}(x)}{dx} \left[ 1 - \frac{E_{imp}}{E_{bone}} \right] \quad (37) $$

with the gradient of the second-order area moment reading as

$$ \begin{aligned} \frac{dI_{z,bone}(x)}{dx} &= -\frac{dI_{z,imp}(x)}{dx} \\ &= \frac{x}{2} \left[ (R_{bone})^4 - (R_{bone} - t_{bone})^4 \right] \times \\ & \quad \times \left[ \left( R_{bone} - \frac{t_{bone}}{2} \right) \sqrt{1 - \frac{(R_{imp})^2 - x^2}{(R_{bone} - \frac{t_{bone}}{2})^2}} \right]^{-1} \times \\ & \quad \times \left( \frac{2 - \frac{(R_{imp})^2 - x^2}{(R_{bone} - \frac{t_{bone}}{2})^2}}{\sqrt{(R_{imp})^2 - x^2}} - \frac{\sqrt{(R_{imp})^2 - x^2}}{(R_{bone} - \frac{t_{bone}}{2})^2} \right) \end{aligned} \quad (38) $$

Finally, in Eq. (35), the normal stress gradients arising from a change in the elastic modulus are given by the following mathematical expression

$$ \left( \frac{\partial \sigma_{xx}^{M_z}}{\partial x} \right)^{dE} = - \frac{\partial E(x, s)}{\partial x} \frac{1}{E_{bone}} \frac{y(s) M_z(x)}{I_{z,bone}(x) + \frac{E_{imp}}{E_{bone}} I_{z,imp}(x)} \quad (39) $$

with the gradient of the mechanical properties as shown in Eq. (33), together with Eq. (34).

### 2.3.4. Axial gradients of normal stresses due to axial loading

As regards the gradients of axial normal stress due to axial force, Eq. (5) implies that

$$ \left( \frac{\partial \sigma_{xx}^{N}}{\partial x} \right) = \left( \frac{\partial \sigma_{xx}^{N}}{\partial x} \right)^{dA} + \left( \frac{\partial \sigma_{xx}^{N}}{\partial x} \right)^{dE} \quad (40) $$

with the new expression analogous to those introduced below Eqs. (28) and (35). In Eq. (40), the normal stress gradients arising from changes in the cross-sectional
geometrical properties are given by the following mathematical expression

$$ \left(\frac{\partial \sigma_{xx}^N}{\partial x}\right)^{dA} = - \frac{E(x, s)}{E_{bone}} \frac{N(x)}{[A_{bone}(x) + A_{imp}(x)(E_{imp}/E_{bone})]^2} \times \frac{dA_{bone}(x)}{dx} \left[1 - \frac{E_{imp}}{E_{bone}}\right] \quad (41) $$

with the gradient of the area with x reading as

$$ \frac{dA_{bone}(x)}{dx} = -\frac{dA_{imp}(x)}{dx} = \frac{4 t_{bone} x}{\sqrt{(R_{imp})^2 - x^2} \sqrt{1 - \frac{(R_{imp})^2 - x^2}{(R_{bone} - \frac{t_{bone}}{2})^2}}} \quad (42) $$

Moreover, in Eq. (40), the normal stress gradients arising from a change in the elastic modulus are given by the following mathematical expression

$$ \left(\frac{\partial \sigma_{xx}^N}{\partial x}\right)^{dE} = \frac{\partial E(x, s)}{\partial x} \frac{1}{E_{bone}} \frac{N(x)}{A_{bone}(x) + A_{imp}(x)(E_{imp}/E_{bone})} \quad (43) $$

with the gradient of the mechanical properties as shown in Eq. (33), together with Eq. (34).

## 2.4. Shear stresses from circumferential integration of normal stress gradients

### 2.4.1. Local equilibrium conditions fulfilled by different stress portions

The local equilibrium conditions according to (25) are not only required to hold for the entire stress state prevailing at a particular location in the bone-implant structure, but also for individual stress portions making up the total stress, according to Eqs. (28), (35), and (40). As regards the stress portions associated with bending around the y-axis and changes in the corresponding bending moment, local equilibrium according to Eq. (25) implies the following circumferential gradients of shear stresses

$$ \left(\frac{\partial \sigma_{xs}^{M_y}}{\partial s}\right)^{dM_y} = -\left(\frac{\partial \sigma_{xx}^{M_y}}{\partial x}\right)^{dM_y} \quad (44) $$

As regards the stress portions associated with bending around the z-axis and changes in the corresponding bending moment, local equilibrium according to Eq. (25) implies the following circumferential gradients of shear stresses

$$ \left(\frac{\partial \sigma_{xs}^{M_z}}{\partial s}\right)^{dM_z} = -\left(\frac{\partial \sigma_{xx}^{M_z}}{\partial x}\right)^{dM_z} \quad (45) $$

As regards the stress portions associated with geometrical and material property changes under bending and axial loading, respectively, local equilibrium according
to Eq. (25) implies the following circumferential gradients of shear stresses

$$ \left(\frac{\partial \sigma_{xs}^{M_y}}{\partial s}\right)^{dA+dE} = - \left(\frac{\partial \sigma_{xx}^{M_y}}{\partial x}\right)^{dA} - \left(\frac{\partial \sigma_{xx}^{M_y}}{\partial x}\right)^{dE} \quad (46) $$

$$ \left(\frac{\partial \sigma_{xs}^{M_z}}{\partial s}\right)^{dA+dE} = - \left(\frac{\partial \sigma_{xx}^{M_z}}{\partial x}\right)^{dA} - \left(\frac{\partial \sigma_{xx}^{M_z}}{\partial x}\right)^{dE} \quad (47) $$

$$ \left(\frac{\partial \sigma_{xs}^{N}}{\partial s}\right)^{dA+dE} = - \left(\frac{\partial \sigma_{xx}^{N}}{\partial x}\right)^{dA} - \left(\frac{\partial \sigma_{xx}^{N}}{\partial x}\right)^{dE} \quad (48) $$

### 2.4.2. Shear stress distribution arising from an axial gradient of the bending moment around the y-axis

According to Eqs. (44), (29), (23) and (10), the circumferential gradients of shear stresses arising from an axial gradient of the bending moment around the y-axis are given by the following mathematical expression

$$ \left(\frac{\partial \sigma_{xs}^{M_y}}{\partial s}\right)^{dM_y} = - \frac{S_z}{\mathcal{I}_y} \left(R_{bone} - \frac{t_{bone}}{2}\right) \times \sin\left[\frac{s}{R_{bone} - \frac{t_{bone}}{2}}\right] \times \frac{1}{E_{bone}} \left[E_{imp} - (E_{imp} - E_{bone}) \times \{H(\xi_1) - H(\xi_2) + H(\xi_3) - H(\xi_4)\}\right] \quad (49) $$

where we introduced the following abbreviations for the arguments of the Heaviside function

$$ \xi_1 = s - \gamma_{int}(x) \times \left(R_{bone} - \frac{t_{bone}}{2}\right) \quad (50) $$

$$ \xi_2 = s - [\pi - \gamma_{int}(x)] \times \left(R_{bone} - \frac{t_{bone}}{2}\right) \quad (51) $$

$$ \xi_3 = s - [\pi + \gamma_{int}(x)] \times \left(R_{bone} - \frac{t_{bone}}{2}\right) \quad (52) $$

$$ \xi_4 = 2\pi - \gamma_{int}(x) \times \left(R_{bone} - \frac{t_{bone}}{2}\right) \quad (53) $$

as well as a “composite second-order area moment around the y-axis” of the form

$$ \mathcal{I}_y(x) = I_{y,bone}(x) + I_{y,imp}(x) \frac{E_{imp}}{E_{bone}} \quad (54) $$

Integration of the shear stress gradients of Eq.(49), together with Eqs. (49) – (54), over the circumferential line of the thin-walled bone profile and its extensions into the implant portion of the cross section, yields the shear stresses arising from the axial
gradient of the bending moment around the y-axis, in the following form,

$$ (\sigma_{xs}^{M_y})^{dM_y}(x, s) = \int_{0}^{s} \left( \frac{\partial \sigma_{xs}^{M_y}}{\partial \tilde{s}} \right)^{dM_y} (x, \tilde{s}) d\tilde{s} + (C^{M_y})^{dM_y}(x) \quad (55) $$

Thereby, the integral expression occurring in Eq. (55) reads as

$$ \left[ \int_{0}^{s} \left( \frac{\partial \sigma_{xs}^{M_y}}{\partial \tilde{s}} \right)^{dM_y} d\tilde{s} \right] (x, s) = \frac{-S_z}{\mathcal{I}_y(x)} \times \frac{1}{E_{bone}} \left( R_{bone} - \frac{t_{bone}}{2} \right)^2 \times \\ \times \left[ E_{imp} (1 - \cos[\gamma(s)]) + (E_{imp} - E_{bone}) \times \right. \\ \times \left\{ (\cos[\gamma(s)] - \cos[\gamma_{int}(x)]) \times [H(\xi_1) - H(\xi_4)] \times \right. \\ \times \left. (\cos[\gamma(s)] + \cos[\gamma_{int}(x)]) \times [-H(\xi_2) + H(\xi_3)] \right\} \right] \quad (56) $$

with $\gamma(s)$ according to Eq. (8), $\gamma_{int}(x)$ according to Eq. (11), $\xi_1$ to $\xi_4$ according to Eqs. (50)–(53), and $\mathcal{I}_y$ according to Eq. (54). The integration constant $(C)^{dM_y}(x)$ occurring in Eq. (55) follows from the condition that shear stresses arising from an axial gradient of the bending moment around the y-axis do not result in a torsional moment (around the x-axis), i.e. from

$$ (M_x^{M_y})^{dM_y} = \int_A -z (\sigma_{xy}^{M_y})^{dM_y} + y(\sigma_{xz}^{M_y})^{dM_y} dA \\ \approx t_{bone} \left( R_{bone} - \frac{t_{bone}}{2} \right) \int_{s=0}^{(R_{bone}-t_{bone}/2)} (\sigma_{xs}^{M_y})^{dM_y} ds = 0 \quad (57) $$

yielding

$$ (C^{M_y})^{dM_y}(x) = \frac{S_z (R_{bone} - \frac{t_{bone}}{2})^2}{\mathcal{I}_y(x) E_{bone}} \left[ E_{imp} - (E_{imp} - E_{bone}) \cos[\gamma_{int}(x)] \right] \quad (58) $$

The shear stresses according to Eqs. (55), (56), and (58) also fulfill two conditions in line with Eq. (3), namely that shear stresses arising from bending gradients around the y-axis do not result in shear forces in the y-direction,

$$ (S_y^{M_y})^{dM_y} = \int_A (\sigma_{xy}^{M_y})^{dM_y} dA \\ \approx -t_{bone} \int_{s=0}^{2\pi(R_{bone}-t_{bone}/2)} \sin[\gamma(s)] \times (\sigma_{xs}^{M_y})^{dM_y} ds = 0 \quad (59) $$
and that the integration, over the cross section, of shear stresses arising from a gradient in the bending moment around the y-axis, yields the shear force in z-direction, i.e.

$$ (S_z^{M_y})^{dM_y} = \int_A (\sigma_{xz}^{M_y})^{dM_y} dA \\ \approx t_{bone} \int_{s=0}^{2\pi(R_{bone}-t_{bone}/2)} \cos[\gamma(s)] \times (\sigma_{xs}^{M_y})^{dM_y} ds \approx S_z \quad (60) $$

### 2.4.3. Shear stress distribution arising from an axial gradient of the bending moment around the z-axis

According to Eqs. (45), (36), (23) and (9), the circumferential gradients of shear stresses arising from an axial gradient of the bending moment around the z-axis are given in the following mathematical form

$$ \left(\frac{\partial \sigma_{xs}^{M_z}}{\partial s}\right)^{dM_z} = \frac{-S_y}{\mathcal{I}_z(x)} \left(R_{bone} - \frac{t_{bone}}{2}\right) \times \cos\left[\frac{s}{R_{bone} - \frac{t_{bone}}{2}}\right] \times \\ \times \frac{1}{E_{bone}} \left[E_{imp} - (E_{imp} - E_{bone}) \times \right. \\ \times \left\{H(\xi_1) - H(\xi_2) + H(\xi_3) - H(\xi_4)\right\} \quad (61) $$

where the abbreviations for the arguments of the Heaviside function still follow Eqs. (50)–(53) and where we introduced a “composite second-order area moment around the z-axis” of the form

$$ \mathcal{I}_z(x) = I_{z,bone}(x) + I_{z,imp}(x) \frac{E_{imp}}{E_{bone}} \quad (62) $$

Integration of the shear stress gradients according to Eq. (61), together with Eq. (62), over the circumferential line of the thin-walled bone profile and its extensions into the implant portion of the cross section, yields the shear stresses arising from the axial gradient of the bending moment around z, in the following form,

$$ (\sigma_{xs}^{M_z})^{dM_z}(x, s) = \int_0^s \left(\frac{\partial \sigma_{xs}^{M_z}}{\partial \tilde{s}}\right)^{dM_z}(x, \tilde{s}) d\tilde{s} + (C^{M_z})^{dM_z}(x) \quad (63) $$
Thereby, the integral expression occurring in Eq. (63) reads as

$$ \left[ \int_{0}^{s} \left( \frac{\partial \sigma_{xs}^{M_z}}{\partial \tilde{s}} \right)^{dM_z} d\tilde{s} \right] (x, s) = \frac{-S_y}{\mathcal{I}_z(x)} \frac{1}{E_{bone}} \left( R_{bone} - \frac{t_{bone}}{2} \right)^2 \times \\ \times \left[ E_{imp} \sin[\gamma(s)] + (E_{imp} - E_{bone}) \times \right. \\ \times \left\{ (\sin[\gamma(s)] - \sin[\gamma_{int}(x)]) \times [-H(\xi_1) + H(\xi_2)] \times \right. \\ \times \left. (\sin[\gamma(s)] + \sin[\gamma_{int}(x)]) \times [-H(\xi_3) + H(\xi_4)] \right\} \right] \quad (64) $$

with $\gamma(s)$ according to (8), $\gamma_{int}(x)$ according to (11), $\xi_1$ to $\xi_4$ according to (50)–(53), and $\mathcal{I}_y$ according to (54). The integration constant $(C^{M_z})^{dM_z}(x)$ follows from the condition that shear stresses arising from bending changes around the z-axis do not result in a torsional moment (around the x-axis), i.e. from

$$ (M_x^{M_z})^{dM_z} = \int_A^{-z} (\sigma_{xy}^{M_z})^{dM_z} + y(\sigma_{xz}^{M_z})^{dM_z} dA \\ \approx t_{bone} \left(R_{bone} - \frac{t_{bone}}{2}\right) \int_0^{(R_{bone} - t_{bone})/2} (\sigma_{xs}^{M_z})^{dM_z} ds = 0 \quad (65) $$

yielding

$$ (C^{M_z})^{dM_z}(x) = 0 \quad (66) $$

The shear stresses according to Eq. (63), together with Eq. (64) and (66), also fulfill two conditions in line with Eq. (3), namely that shear stresses arising from an axial gradient of the bending moment around the z-axis do not result in shear forces in the y-direction,

$$ (S_z^{M_z})^{dM_z} = \int_A (\sigma_{xz}^{M_z})^{dM_z} dA \\ \approx -t_{bone} \int_0^{2\pi(R_{bone}-t_{bone}/2)} \sin[\gamma(s)] (\sigma_{xs}^{M_z})^{dM_z} ds = 0 \quad (67) $$

and that the integration, over the cross section, of the shear stresses arising from an axial gradient of the bending moment around the z-axis yields the shear force in the y-direction, i.e.

$$ (S_y^{M_z})^{dM_z} = \int_A (\sigma_{xy}^{M_z})^{dM_z} dA \\ \approx t_{bone} \int_0^{2\pi(R_{bone}-t_{bone}/2)} \cos[\gamma(s)] (\sigma_{xs}^{M_z})^{dM_z} ds = S_y \quad (68) $$
### 2.4.4. Shear stress distribution arising from axial gradients of cross-sectional properties under bending around the y-axis

According to Eqs. (46), (30), (32), and (10), the circumferential gradients of shear stresses arising from axial gradients of the cross-sectional geometrical properties and of the mechanical properties of a bone-implant compound cross section loaded by a bending moment around the y-axis are given by the following mathematical expression

$$ \left(\frac{\partial \sigma_{xs}^{M_y}}{\partial s}\right)^{dA+dE} (x, s) = \frac{1}{E_{bone}} \frac{M_y(x)z(s)}{\mathcal{I}_y(x)} \times \\ \times \left[ \frac{1}{\mathcal{I}_y(x)} \frac{dI_{y,bone}(x)}{dx} \left[ 1 - \frac{E_{imp}}{E_{bone}} \right] \times E(x, s) - \frac{\partial E(x, s)}{\partial x} \right] = \\ = \frac{1}{E_{bone}} \frac{M_y(x)}{\mathcal{I}_y(x)} \left( R_{bone} - \frac{t_{bone}}{2} \right) \sin \left[ \frac{s}{R_{bone} - \frac{t_{bone}}{2}} \right] \times \\ \times \left[ \frac{1}{\mathcal{I}_y(x)} \frac{dI_{y,bone}(x)}{dx} \left[ 1 - \frac{E_{imp}}{E_{bone}} \right] \left[ E_{imp} - (E_{imp} - E_{bone}) \right] \times \right. \\ \left. \times \left\{ H(\xi_1) - H(\xi_2) + H(\xi_3) - H(\xi_4) \right\} \right] \\ - \frac{d\gamma_{int}(x)}{dx} \left( R_{bone} - \frac{t_{bone}}{2} \right) (E_{imp} - E_{bone}) \times \\ \times \left\{ \delta(\xi_1) + \delta(\xi_2) + \delta(\xi_3) + \delta(\xi_4) \right\} \right] \quad (69) $$

with the composite second-order area moment $\mathcal{I}_y$ being given by Eq. (54). Integration of the shear gradients of Eq. (69), over the circumferential line of the bone profile and its extensions into the implant portions of the cross section, yields the shear stresses arising from the axial gradients of the cross-sectional geometrical properties and of the mechanical properties of a bone-implant compound cross-section loaded by a bending moment around the y-axis as

$$ (\sigma_{xs}^{M_y})^{dA+dE} (x, s) = (C^{M_y})^{dA+dE} (x) + \int_{0}^{s} \left( \frac{\partial \sigma_{xs}^{M_y}}{\partial \tilde{s}} \right)^{dA+dE} (x, \tilde{s}) d\tilde{s} \quad (70) $$
Thereby, the integral expression in Eq. (70) reads as

$$ \left[ \int_{0}^{s} \left( \frac{\partial \sigma_{xs}^{M_y}}{\partial \tilde{s}} \right)^{dA+dE} d\tilde{s} \right] (x, s) = \frac{M_y}{\mathcal{I}_y(x)} \frac{1}{E_{bone}} \left( R_{bone} - \frac{t_{bone}}{2} \right)^2 \times \\ \times \left( \frac{1}{\mathcal{I}_y(x)} \frac{dI_{y,bone}(x)}{dx} \left[ 1 - \frac{E_{imp}}{E_{bone}} \right] \times \right. \\ \times \left[ -E_{imp} \cos[\gamma(s)] + (E_{imp} - E_{bone}) \times \right. \\ \times \left\{ (\cos[\gamma(s)] - \cos[\gamma_{int}(x)]) [H(\xi_1) - H(\xi_4)] \right. \\ + (\cos[\gamma(s)] + \cos[\gamma_{int}(x)]) [H(\xi_2) - H(\xi_3)] \right\} \right] \\ + (E_{imp} - E_{bone}) \sin[\gamma_{int}(x)] \frac{d\gamma_{int}(x)}{dx} \times \\ \times \left\{ -H(\xi_1) - H(\xi_2) + H(\xi_3) + H(\xi_4) \right\} \quad (71) $$

with $\gamma(s)$ according to (8), $\gamma_{int}(x)$ according to (11), $\xi_1$ to $\xi_4$ according to (50)–(53), and $\mathcal{I}_y$ according to (54). The integration constant $(C^{M_y})^{dA+dE}$ ($x$) follows from the condition that shear stresses arising from axial gradients of the cross-sectional geometrical properties and of the mechanical properties of a bone-implant compound cross section loaded by a bending moment around the z-axis do neither result in a torsional moment (around the x-axis), i.e. from

$$ (M_x^{M_y})^{dA+dE} = \int_A -z (\sigma_{xy}^{M_y})^{dA+dE} + y (\sigma_{xz}^{M_y})^{dA+dE} dA \\ \approx t_{bone} \left(R_{bone} - \frac{t_{bone}}{2}\right) \int_0^{2\pi(R_{bone}-t_{bone}/2)} (\sigma_{xs}^{M_y})^{dA+dE} dA = 0 \quad (72) $$

nor in shear forces in y- and z-direction

$$ (S_y^{M_y})^{dA+dE} = \int_A (\sigma_{xy}^{M_y})^{dA+dE} dA \\ \approx t_{bone} \int_0^{2\pi(R_{bone}-t_{bone}/2)} \cos[\gamma(s)] (\sigma_{xs}^{M_y})^{dA+dE} ds = 0 \quad (73) $$

$$ (S_z^{M_y})^{dA+dE} = \int_A (\sigma_{xz}^{M_y})^{dA+dE} dA \\ \approx t_{bone} \int_0^{2\pi(R_{bone}-t_{bone}/2)} -\sin[\gamma(s)] (\sigma_{xs}^{M_y})^{dA+dE} ds = 0 \quad (74) $$
yielding

$$ (C^{My})^{dA+dE}(x) = \frac{M_y}{\mathcal{I}_y(x)} \frac{1}{E_{bone}} \left(R_{bone} - \frac{t_{bone}}{2}\right)^2 \left[1 - \frac{E_{imp}}{E_{bone}}\right] \times \left[ \frac{E_{imp} - E_{bone}}{\mathcal{I}_y(x)} \frac{dI_{y,bone}(x)}{dx} \cos[\gamma_{int}(x)] - E_{bone} \frac{d\gamma_{int}(x)}{dx} \sin[\gamma_{int}(x)] \right] \quad (75) $$

### 2.4.5. Shear stress distribution arising from axial gradients of cross-sectional properties under bending around the z-axis

According to Eqs. (47), (37), (39) and (9), the circumferential gradients of shear stresses arising from axial gradients of the cross-sectional geometrical properties and of mechanical properties of a bone-implant compound cross section loaded by a bending moment around the z-axis are given by the following mathematical expression

$$ \left(\frac{\partial \sigma_{xs}^{M_z}}{\partial s}\right)^{dA+dE}(x, s) = -\frac{1}{E_{bone}} \frac{M_z(x) y(s)}{\mathcal{I}_z(x)} \times \\ \times \left[ \frac{1}{\mathcal{I}_z(x)} \frac{dI_{z,bone}(x)}{dx} \left[ 1 - \frac{E_{imp}}{E_{bone}} \right] E(x, s) - \frac{\partial E(x, s)}{\partial x} \right] \\ = -\frac{M_z(x)}{\mathcal{I}_z(x)} \frac{1}{E_{bone}} \left(R_{bone} - \frac{t_{bone}}{2}\right) \cos\left[\frac{s}{R_{bone} - \frac{t_{bone}}{2}}\right] \times \\ \times \left[ \frac{1}{\mathcal{I}_z(x)} \frac{dI_{z,bone}(x)}{dx} \left[ 1 - \frac{E_{imp}}{E_{bone}} \right] \left[ E_{imp} - (E_{imp} - E_{bone}) \times \right. \right. \\ \left. \left. \times \{H_1(x, s) - H_2(x, s) + H_3(x, s) - H_4(x, s)\} \right] \right. \\ \left. - \frac{d\gamma_{int}(x)}{dx} \left(R_{bone} - \frac{t_{bone}}{2}\right) (E_{imp} - E_{bone}) \times \right. \\ \left. \times \{\delta_1(x, s) + \delta_2(x, s) + \delta_3(x, s) + \delta_4(x, s)\} \right] \quad (76) $$

with the composite second-order area moment $\mathcal{I}_z$ being given by Eq. (62). Integration of the shear stress gradients according to Eq. (76), over the circumferential line of the bone profile and its extensions into the implant portion of the cross section, yields the shear stresses arising from the axial gradients of the cross-sectional geometrical properties and of the mechanical properties of a bone-implant compound cross-section loaded by a bending moment around the z-axis as

$$ (\sigma_{xs}^{M_z})^{dA+dE}(x, s) = (C^{Mz})^{dA+dE}(x) + \int_{0}^{s} \left(\frac{\partial \sigma_{xs}^{M_z}}{\partial \tilde{s}}\right)^{dA+dE}(x, \tilde{s}) d\tilde{s} \quad (77) $$
Thereby, the integral expression in Eq. (77) reads as

$$ \left[ \int_{0}^{s} \left( \frac{\partial \sigma_{xs}^{M_z}}{\partial \tilde{s}} \right)^{dA+dE} d\tilde{s} \right] (x, s) = \frac{M_z(x)}{\mathcal{I}_z(x)} \frac{1}{E_{bone}} \left( R_{bone} - \frac{t_{bone}}{2} \right)^2 \times \\ \times \left( \frac{1}{\mathcal{I}_z(x)} \frac{dI_{z,bone}(x)}{dx} \left[ 1 - \frac{E_{imp}}{E_{bone}} \right] \times \right. \\ \times \left[ -E_{imp} \sin[\gamma(s)] + (E_{imp} - E_{bone}) \times \right. \\ \times \left\{ (\sin[\gamma(s)] - \sin[\gamma_{int}(x)]) [H(\xi_1) - H(\xi_2)] \right. \\ \times \left. + (\sin[\gamma(s)] + \sin[\gamma_{int}(x)]) [H(\xi_3) - H(\xi_4)] \right\} \right] \\ + (E_{imp} - E_{bone}) \cos[\gamma_{int}(x)] \frac{d\gamma_{int}(x)}{dx} \times \\ \times \left[ H(\xi_1) - H(\xi_2) - H(\xi_3) + H(\xi_4) \right] \tag{78} $$

with $\gamma(s)$ according to in (8), $\gamma_{int}(x)$ according to in (11), $\xi_1$ to $\xi_4$ according to (50)–(53), and $\mathcal{I}_y$ according to (54). The integration constant $(C^{Mz})^{dA+dE}$ ($x$) follows from the condition that shear stresses arising from axial gradients in the cross-sectional geometrical properties and of the mechanical properties of a bone-implant compound cross-section loaded by a bending moment around the $z$-axis do neither result in a torsional moment (around the $x$-axis), i.e. from

$$ \begin{aligned} (M_x^{M_z})^{dA+dE} &= \int_A -z (\sigma_{xy}^{M_z})^{dA+dE} + y(\sigma_{xz}^{M_z})^{dA+dE} dA \\ &\approx t_{bone} \left(R_{bone} - \frac{t_{bone}}{2}\right) \int_0^{2\pi(R_{bone}-t_{bone}/2)} (\sigma_{xs}^{M_z})^{dA+dE} ds = 0 \end{aligned} \tag{79} $$

nor in shear forces in $y$- and $z$-direction

$$ \begin{aligned} (S_y^{M_z})^{dA+dE} &= \int_A (\sigma_{xy}^{M_z})^{dA+dE} dA \\ &\approx t_{bone} \int_A \cos[\gamma(s)] (\sigma_{xs}^{M_z})^{dA+dE} ds = 0 \end{aligned} \tag{80} $$

$$ \begin{aligned} (S_z^{M_z})^{dA+dE} &= \int_A (\sigma_{xz}^{M_z})^{dA+dE} dA \\ &\approx t_{bone} \int_0^{2\pi(R_{bone}-t_{bone}/2)} -\sin[\gamma(s)] (\sigma_{xs}^{M_z})^{dA+dE} ds = 0 \end{aligned} \tag{81} $$

yielding

$$ (C^{Mz})^{dA+dE} (x) = 0 \tag{82} $$
### 2.4.6. Shear stress distribution arising from axial gradients of cross-sectional properties under axial loading

According to Eqs. (48), (41) and (43), the circumferential gradients of shear stresses arising from axial gradients of the cross-sectional geometrical properties and of the mechanical properties of a bone-implant compound cross section loaded by a constant axial force are given by the following mathematical expression

$$ \left(\frac{\partial \sigma_{xs}^{N}}{\partial s}\right)_{(x, s)} = \frac{1}{E_{bone}} \frac{N(x)}{A(x)} \times \\ \times \left( \frac{1}{A(x)} \frac{dA_{bone}(x)}{dx} \left[ 1 - \frac{E_{imp}}{E_{bone}} \right] \left[ E_{imp} - (E_{imp} - E_{bone}) \times \right. \right. \\ \left. \left. \times \{H(\xi_1) - H(\xi_2) + H(\xi_3) - H_4(\xi_4)\} \right] \right. \\ \left. \times - \frac{d\gamma_{int}(x)}{dx} \left(R_{bone} - \frac{t_{bone}}{2}\right) (E_{imp} - E_{bone}) \times \right. \\ \left. \times \{\delta(\xi_1) + \delta(\xi_2) + \delta(\xi_3) + \delta(\xi_4)\} \right) \quad (83) $$

where we introduced a “composite area” of the form

$$ A(x) = A_{bone}(x) + A_{imp}(x) \frac{E_{imp}}{E_{bone}} \quad (84) $$

Integration of the shear gradients according to Eq. (83), over the circumferential line of the bone profile and its extensions into the implant portion of the cross section, yields the shear stresses arising from the changes in the cross-sectional geometrical properties and of mechanical properties of a bone-implant compound cross-section loaded by a constant axial force as

$$ \left(\sigma_{xs}^{N}\right)^{dA+dE}(x, s) = \left(C^{N}\right)^{dA+dE}(x) + \int_{0}^{s} \left(\frac{\partial \sigma_{xs}^{N}}{\partial \tilde{s}}\right)^{dA+dE}(x, \tilde{s}) d\tilde{s} \quad (85) $$
Thereby, the integral expression occurring in Eq. (85) reads as

$$ \left[ \int_{0}^{s} \left( \frac{\partial \sigma_{xs}^{N}}{\partial \tilde{s}} \right)^{dA+dE} d\tilde{s} \right] (x, s) = \frac{N(x)}{\mathcal{A}(x)} \frac{1}{E_{bone}} \times \\ \times \left( \frac{1}{\mathcal{A}(x)} \frac{dA_{bone}(x)}{dx} \left[ 1 - \frac{E_{imp}}{E_{bone}} \right] \times \right. \\ \times \left[ E_{imp} \times s + (E_{imp} - E_{bone}) \times \right. \\ \times \left\{ -\xi_{1}H(\xi_{1}) + \xi_{2}H(\xi_{2}) - \xi_{3}H(\xi_{3}) + \xi_{4}H(\xi_{4}) \right\} \\ - \frac{d\gamma_{int}(x)}{dx} \left( R_{bone} - \frac{t_{bone}}{2} \right) (E_{imp} - E_{bone}) \times \\ \times \left\{ H(\xi_{1}) + H(\xi_{2}) + H(\xi_{3}) + H(\xi_{4}) \right\} \right] \quad (86) $$

with $\gamma(s)$ according to (8), $\gamma_{int}(x)$ according to (11), $\xi_1$ to $\xi_4$ according to (50)–(53), and $\mathcal{A}$ as in (84). The integration constant $(C^N)^{dA+dE}(x)$ follows from the condition that shear stresses arising from axial gradients in the cross-sectional geometrical properties and of the mechanical properties of a bone-implant compound cross-section loaded by a bending moment around the z-axis do neither result in a torsional moment (around the x-axis), i.e. from

$$ (M_x^N)^{dA+dE} = \int_A -z (\sigma_{xy}^N)^{dA+dE} + y (\sigma_{xz}^N)^{dA+dE} dA \\ \approx t_{bone} \left(R_{bone} - \frac{t_{bone}}{2}\right) \int_0^{2\pi(R_{bone}-t_{bone}/2)} (\sigma_{xs}^N)^{dA+dE} ds = 0 \quad (87) $$

nor in shear forces in y- and z-direction

$$ (S_y^N)^{dA+dE} = \int_A (\sigma_{xy}^N)^{dA+dE} dA \\ \approx t_{bone} \int_0^{2\pi(R_{bone}-t_{bone}/2)} \cos[\gamma(s)] (\sigma_{xs}^N)^{dA+dE} ds = 0 \quad (88) $$

$$ (S_z)^{N,dA+dE} = \int_A (\sigma_{xz}^N)^{dA+dE} dA \\ \approx \int_0^{2\pi(R_{bone}-t_{bone}/2)} \left( -\sin[\gamma(s)](\sigma_{xs}^N)^{dA+dE} \right) dA = 0 \quad (89) $$

yielding

$$ (C^N)^{dA+dE}(x) = \frac{-N}{\mathcal{A}(x)} \frac{1}{E_{bone}} \times \left[ \frac{E_{bone} \pi + 2(E_{imp} - E_{bone})\gamma_{int}(x)}{\mathcal{A}(x)} \times \frac{dA_{bone}(x)}{dx} \left[ 1 - \frac{E_{imp}}{E_{bone}} \right] - 2(E_{imp} - E_{bone}) \frac{d\gamma_{int}(x)}{dx} \right] = 0 \quad (90) $$
## 2.5. Shear stresses arising from torsional loading

As regards shear stresses arising from torsional loading, we adopt the classical modeling approach. Accordingly, the local equilibrium condition for shear stresses $\sigma_{xs}^{M_x}$ arising from a torsional moment $M_x$, reads as

$$ \frac{\partial \sigma_{xs}^{M_x}}{\partial s} = 0 \rightarrow \sigma_{xs}^{M_x}(s, x = \text{const.}) = \text{const.} \qquad (91) $$

and these shear stresses obey the classical relation ([19], 2011, page 414)

$$ \sigma_{xs}^{M_x} = \frac{M_x (R_{bone} - \frac{t_{bone}}{2})}{I_p} \qquad (92) $$

with the second polar moment of area reading as

$$ I_p = \frac{\pi}{2} [(R_{bone})^4 - (R_{bone} - t_{bone})^4] \qquad (93) $$

## 2.6. Higher-order stress average as measure for implant and bone loading

In order to study the effect of design parameter changes, such as the radius or the elastic modulus of the implant, on the stress levels encountered in the bone and implant portions of the femur-implant compound structure, we resort to the concept of higher-order stress averages [23,24], which have turned out as very useful measures in the context of strength upscaling in micromechanics [25], with applications to bone biomaterials [26], metallic implant materials [27], wood [28,29], and concrete [30–32]. As concerns appropriate measures for the load bearing capacity, the stress deviator $\sigma_d$ is a particularly popular choice,

$$ \sigma_d = \sigma - \frac{1}{3} \text{tr} \sigma \qquad (94) $$
with the corresponding higher-order stress average over the bone and implant portions in a cross section of the organ-implant compound structure being given as [24,28,30,33]

$$ 
\begin{aligned}
\mathfrak{s}_{\text{bone}} &= \sqrt{\frac{1}{A_{\text{bone}}} \int_{A_{\text{bone}}} \frac{1}{2} \boldsymbol{\sigma}_d : \boldsymbol{\sigma}_d dA} \\
&\approx \sqrt{\frac{1}{2\pi - 4\gamma_{\text{int}}} \left( \int_{\gamma_{\text{int}}}^{\pi - \gamma_{\text{int}}} \left[ (\sigma_{xs})^2 + \frac{(\sigma_{xx})^2}{3} \right] d\gamma + \int_{\pi + \gamma_{\text{int}}}^{2\pi - \gamma_{\text{int}}} \left[ (\sigma_{xs})^2 + \frac{(\sigma_{xx})^2}{3} \right] d\gamma \right)} \quad (95)
\end{aligned}
$$

$$ 
\begin{aligned}
\mathfrak{s}_{\text{imp}} &= \sqrt{\frac{1}{A_{\text{imp}}} \int_{A_{\text{imp}}} \frac{1}{2} \boldsymbol{\sigma}_d : \boldsymbol{\sigma}_d dA} \\
&\approx \sqrt{\frac{1}{4\gamma_{\text{int}}} \left( \int_{0}^{\gamma_{\text{int}}} \left[ (\sigma_{xs})^2 + \frac{(\sigma_{xx})^2}{3} \right] d\gamma + \int_{\pi - \gamma_{\text{int}}}^{\pi + \gamma_{\text{int}}} \left[ (\sigma_{xs})^2 + \frac{(\sigma_{xx})^2}{3} \right] d\gamma \right.} \\
&\qquad \left. + \int_{2\pi - \gamma_{\text{int}}}^{2\pi} \left[ (\sigma_{xs})^2 + \frac{(\sigma_{xx})^2}{3} \right] d\gamma \right)} \quad (96)
\end{aligned}
$$

with the stress tensor following Eq. (1), together with Eqs. (26) and (27), whereby the normal and shear stress components result from summing up contributions from normal forces, bending and torsional moments, and those from changes in bending moments and cross-sectional as well as material properties

$$ \sigma_{xx} = \sigma_{xx}^{N} + \sigma_{xx}^{M_y} + \sigma_{xx}^{M_z} \quad (97) $$

$$ \sigma_{xs} = (\sigma_{xs}^{N})^{dA+dE} + (\sigma_{xs}^{M_y})^{dA+dE} + (\sigma_{xs}^{M_z})^{dA+dE} + (\sigma_{xs}^{M_y})^{dM_y} + (\sigma_{xs}^{M_z})^{dM_z} \quad (98) $$

### 3. Illustration of geometrical properties and characteristic stress states

While the second-order area moment around the y-axis, of the bone portions of the composite cross sections is hardly influenced by the presence of the implant, see Figure 4(c), this presence reduces the cross-sectional bone area by up to 25%, see Figure 4(c), and the second-order bone area moment around the z-axis by around 50%, see Figure 4(e). According to Eqs. (5) to (7), we therefore expect only slight to moderate effects of the presence of the implant on the axial normal stresses, as compared to the situation without implant, illustrated in Figure 5. At the same time, the axial gradients of all area and area moment measures exhibit singularities at the positions $x = \pm R_{\text{imp}}$, see Figure 4(b,d,f). According to Eqs. (5) to (7), together with Eq. (25), this lets us expect shear stresses which tend towards infinity at these very locations, to be discussed in more detail in the next paragraph.

In quantitative and illustrative terms, the stress distributions across the femur-implant structure can be described as follows: The normal stresses arising from axial loading are scaled by the ratio of the elastic moduli of bone and implant, respectively, irrespective of the thickness of the implant portion in the composite cross sections, see the stress distributions in the third column of Figure 6, as an illustration of the evaluation of Eq. (5) for the loading defined by Eq.(2). A similar scaling is observed for the normal stresses due to bending moments around y- and z-axis, respectively.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Typical trends of cross-sectional area A and area moments *I*<sub>y</sub> and *I*<sub>z</sub> along bone shaft axis around the magnesium implant; i.e. for the geometrical conditions encountered in [11], see Table 1</figcaption></figure>
These bending moments, however, do not only scale with the modulus ratio, but also
with the distance from the geometrical center of the composite cross section, see stress
distributions in the fourth and fifth column of Figure 6, as illustrations of Eqs. (6)
and (7) for the loading defined by Eq.(2). In quantitative terms, the loading given
through Eq.(2) results in normal stresses due to bending around the y-axis, i.e. the
axis parallel to the implant axis, which dominate in the bone portion of the compound
structure, while the normal stresses due to (compressive) axial normal force dominate
in the implant portion. By comparison, the normal stresses due to bending around the
z-axis are significantly smaller. However, gradients of the bending moment around the
z-axis along the long axis of the bone shaft, i.e along the x-axis, result in significant
axial-circumferential shear stresses, see the stress distributions in the fifth column of
Figure 7, as an illustration of Eq. (63), together with Eqs. (64) and (66), for the
loading defined by Eq.(2). These shear stresses are much larger in magnitude than
those arising from torsional loading, see the stress distributions in the third column
of Figure 7, as an illustration of Eq. (92), together with Eq. (93), for the loading
defined by Eq.(2). The shear stresses due to changing bending moments around the
z axis are also much larger than those arising from the change of bending moments
around the y-axis along the long axis of the bone shaft, see the stress distributions in
the fourth column of Figure 7, as an illustration of Eq. (55), together with Eqs. (56)
and (58), for the loading defined by Eq.(2). Still, along the implant-bone interface in
the periphery of the implant, i.e. where x ≈ ±Rimp, the absolute value of all these
shear stresses is significantly exceeded by the magnitude of the shear stresses arising
from changes along the bone shaft axis, in area, in area moment, as well as in elastic
modulus, of the bone- and implant-specific cross-sectional portions, respectively; see
the first and second column of Figure 8. These shear stresses are mainly triggered by
the axial normal force, see third column of Figure 8, they are also due to bending
around the z-axis, see fifth column of Figure 8, while bending around the y-axis plays
a subordinate role, see fourth column of Figure 8. Being zero at the implant center,
x = 0, these cross-sectional-property-driven shear stresses build up fairly quickly when
following the bone-implant interface at the bone side, along the x-axis, see Figure 9.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Normal and shear stresses arising from the loading given in Eq. (2), for the “classical case” of a hollow cylindrical beam without implant ($R_{imp}=0$)</figcaption></figure>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Axial normal stresses $\sigma_{xx}$ arising from normal force (N) and from bending moments around the y- and z-axis ($M_y$ and $M_z$), in different cross sections through the implant (red) and the adjacent bone portions (blue)</figcaption></figure>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Axial-circumferential shear stresses $\sigma_{xs}$ arising from torsional moments ($M_x$), and from changes in the bending moments ($dM_y$ and $dM_z$), in different cross sections through the implant (red) and the adjacent bone portions (blue)</figcaption></figure>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Axial-circumferential shear stresses arising from changes in the cross-section along the beam axis, while being loaded by axial force and bending moments, in different cross sections through the implant (red) and the adjacent bone portions (blue)</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Axial variation of shear stresses arising from constant axial force together with changes in the cross-section, bending moments together with changes in the cross-section, along the bone-side of the bone-implant interface</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Average second invariant of the stress deviator $s_{bone}$ and $s_{imp}$ within selected cross-sections in dependency of $R_{imp}/(R_{bone} - t_{bone})$ for a fixed ratio of $R_{bone}/t_{bone}$; dash-dotted lines refer to geometrical and mechanical conditions given in Table 1</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Average second invariant of the stress deviator $s_{bone}$ and $s_{imp}$ within selected cross-sections in dependency of $E_{imp}/E_{bone}$ for a fixed ratio of $R_{bone}/t_{bone}$; dash-dotted lines refer to geometrical and mechanical conditions given in Table 1</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Alignment of bone-nanostructure preferred orientation with bone-axis from small-angle x-ray scattering tensor tomography, reproduced and adapted from Figure 5 of [5], under license CC BY 4.0</figcaption></figure>

**4. Discussion**

In this paper, we developed a new analytical mechanics model for a bone-implant structure consisting of a femoral shaft and a cylindrical implant positioned orthogonal to the long-bone axis of the femoral shaft: Normal stresses following from classical composite beam theory were combined with local equilibrilium conditions, so as to derive, based on differentiation and integration steps in space, closed-form analytical expressions for shear stress fields throughout the considered preclinical animal model. The most striking feature which distinguishes the results of this new model from the results of classical beam theory applied to the femoral bone shaft only, are the high shear stresses occuring close to the implant surface generators intersecting the bone shaft axis. If these stress singularities are a realistic reflection of dramatic stress re-distributions due to implantation, then they would – according the observations in mechanobiology – cause significant bone remodeling and transformation events in the aforementioned bone regions adjacent to the implant. And indeed, such transformation events are clearly evidenced by tensor tomograms of the surroundings of magnesium implants in murine femora, see Figure 12. These remodeling processes – namely reorientation of bone fibril, away from the bone axis direction, towards a direction which is orthogonal to the bone axis – occur close to the bone-implant interface regions adjacent to the bone axis, and diminish with increasing distance from this axis. Hence, the remodeling events are closely related to the very high shear stresses proposed by our new model, and in this sense, the experimental evidence underlines the principal relevance of the modeling approach. In more detail, the probable mechanobiological phenomena may be sketched as follows: At first, high stress concentrations result in “overuse” of bone, and in local demineralization, delamination, and disintegration of the extracellular nano-composite, i.e. to bone resorption [34], followed by the process of new bone formation, which starts with the deposition of fibrils in a non-mineralized ground substance [35]. If such a material system is subjected to shear, the fibrils do not only undergo deformations, but also rotational rigid body motions, as was shown by dedicated micromechanical models based on Eshelby problems formulated in rate form, together with thermodynamically consistent hypoelastic material modeling at the phase level [36,37]. These rotated fibrils then undergo the well-known vesicle-driven mineralization process [38–40], which affects both intrafibrillar and extrafibrillar spaces [41,42]. Thereby, however, the extrafibrillar mineral portion comprises
not only the majority of the mineral crystals [43–45], but also appears as key driver of important mechanical properties of bone, such as its (visco-)elasticity [46–48] and strength [49,50].

Finally, our approach is naturally characterized by several limitations which may motivate further studies in the future. This naturally relates to geometrical objects exceeding the simple nature of circular cylinders; but, probably more interestingly, to more advanced beam theory approaches, as may be motivated in the following way: The current differentiation and integration of the normal stresses arising from classical composite beam theory result in partially infinitely large shear stresses. While the occurrence of “very large” shear stresses makes great sense in terms of mechanobiology, as discussed further above, these singularities can obviously reflect quantitatively these very large stresses, i.e. they cannot provide numerical values for these stresses. If such a quantitative analysis is desired, the proposed modeling approach needs further refinement. Actually, the occurrence of very high shear stresses, together with very large strain gradients, may propose, as an interesting future step, the consideration of corresponding shear strains in the context of formulating equilibrium conditions at the level of the overall composite beam. Therefore, the virtual fields as given in Eqs. (A1), would need to undergo extensions in their degrees of freedom, so as to allow for significant virtual shear rate gradients [51,52], which follow the geometrical features of the composite cross section. This may lead to additional developments and applications in the field of shear-compliant composite beam theories, as they have been very successfully developed in the context of lightweight and aerospace engineering [53,54].

## 5. Acknowledgements

The authors gratefully acknowledge the financial support provided by the Austrian Science Fund (FWF), in the framework of the project ‘Multiscale changes in bone due to bio-resorbable implants’ (project number: 14409).

## 6. Disclosure statement

No potential conflict of interest was reported by the authors.

## Appendix A. Bending beam equilibrium derived from Principle of Virtual Power

So as to motivate the equilibrium conditions of Eqs. (3), we here resort to the principle of virtual power (PVP) [51,55–57], one of most versatile fundamentals of mechanics with roots going back to ancient times [58,59]. The PVP provides a firm anchorage of beam theory in the 3D continuum mechanics framework. Accordingly, we adopt an arbitrary virtual rigid body rotational motion of a rigid cross section around the y- and z-axes,

$$ \hat{\mathbf{v}}(x, y, z) = \hat{\omega}_y(x)z\mathbf{e}_x - \hat{\omega}_z(x)y\mathbf{e}_x \quad \forall \hat{\omega}_y(x), \hat{\omega}_z(x) \in \mathbb{R} \quad (A1) $$
with $\hat{\omega}_y$ and $\hat{\omega}_z$ being the components of virtual angular velocity vector $\hat{\omega} = \hat{\omega}_y e_y + \hat{\omega}_z e_z$. When considering a situation where this virtual velocity field does not produce any power of external forces, then the virtual power of the internal forces need to vanish in an equilibrated 3D classical continuum. This reads mathematically as

$$P^{\text{int}} = \int_V \boldsymbol{\sigma} : \nabla^S \hat{\mathbf{v}} dV = 0 \qquad (\text{A2})$$

with $\nabla^S$ standing for the symmetrized nabla (or gradient) operator. Insertion of (A1) into (A2) yields

$$\int_{-l/2}^{l/2} \int_A \sigma_{xx} z dA \frac{d\hat{\omega}_y}{dx} dx + \int_{-l/2}^{l/2} \int_A \sigma_{xz} dA \hat{\omega}_y dx - \int_{-l/2}^{l/2} \int_A \sigma_{xy} y dA \frac{d\hat{\omega}_z}{dx} dx - \int_{-l/2}^{l/2} \int_A \sigma_{xy} dA \hat{\omega}_z dx = 0 \quad \forall \hat{\omega}_y(x), \hat{\omega}_z(x) \in \mathbb{R} \quad (\text{A3})$$

The quantities performing power on angular velocities and their gradients are called shear forces and bending moments, respectively; this yields

$$M_y = \int_A \sigma_{xx} z dA \quad M_z = - \int_A \sigma_{xy} y dA \qquad (\text{A4})$$

$$S_y = \int_A \sigma_{xy} dA \quad S_z = \int_A \sigma_{xy} dA \qquad (\text{A5})$$

and insertion of (3) and (A5) into (A3) yields, after integration by parts, that

$$M_y \hat{\omega}_y|_{-l/2}^{l/2} - \int_{-l/2}^{l/2} \left( \frac{dM_y}{dx} - S_z \right) \hat{\omega}_y dx + M_z \hat{\omega}_z|_{-l/2}^{l/2} - \int_{-l/2}^{l/2} \left( \frac{dM_z}{dx} + S_y \right) \hat{\omega}_z dx = 0 \quad \forall \hat{\omega}_y(x), \hat{\omega}_z(x) \in \mathbb{R} \quad (\text{A6})$$

This readily yields the equilibrium conditions (3).

## References

[1] P Mukherjee, S Roy, D Ghosh, and SK Nandi. Role of animal models in biomedical research: a review. *Laboratory Animal Research*, 38(1):18, 2022.

[2] Sonali K Doke and Shashikant C Dhawale. Alternatives to animal testing: A review. *Saudi Pharmaceutical Journal*, 23(3):223-229, 2015.

[3] Christoph Castellani, Richard A Lindtner, Peter Hausbrandt, Elmar Tschegg, Stefanie E Stanzl-Tschegg, Gerald Zanoni, Stefan Beck, and Annelie-Martina Weinberg. Bone-implant interface strength and osseointegration: Biodegradable magnesium alloy versus standard titanium control. *Acta biomaterialia*, 7(1):432-440, 2011.

[4] Tilman A Grünewald, Harald Rennhofer, Bernhard Hesse, Manfred Burghammer, Stefanie E Stanzl-Tschegg, Marine Cotte, Jörg F Löffler, Annelie-Martina Weinberg, and Helga C Lichtenegger. Magnesium from bioresorbable implants: Distribution and impact on the nano-and mineral structure of bone. *Biomaterials*, 76:250-260, 2016.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

[5] Marianne Liebi, Viviane Lutz-Bueno, Manuel Guizar-Sicairos, Bernd M Schönbauer, Johannes Eichler, Elisabeth Martinelli, Jörg F Löffler, Annelie Weinberg, Helga Lichteneger, and Tilman A Grünewald. 3d nanoscale analysis of bone healing around degrading mg implants evaluated by x-ray scattering tensor tomography. *Acta Biomaterialia*, 134:804-817, 2021.

[6] Alexander Hermann, Arman Shojaei, Dirk Steglich, Daniel Höhe, Berit Zeller-Plumhoff, and Christian J Cyron. Combining peridynamic and finite element simulations to capture the corrosion of degradable bone implants and to predict their residual strength. *International journal of mechanical sciences*, 220:107143, 2022.

[7] Ann-Kathrin Gartzke, Stefan Julmi, Christian Klose, Anja-Christina Waselau, Andrea Meyer-Lindenberg, Hans Jürgen Maier, Silke Besdo, and Peter Wriggers. A simulation model for the degradation of magnesium-based bone implants. *Journal of the mechanical behavior of biomedical materials*, 101:103411, 2020.

[8] Mojtaba Barzegari, Di Mei, Sviatlana V Lamaka, and Liesbet Geris. Computational modeling of degradation process of biodegradable magnesium biomaterials. *Corrosion Science*, 190:109674, 2021.

[9] Stefan Scheiner and Christian Hellmich. Stable pitting corrosion of stainless steel as diffusion-controlled dissolution process with a sharp moving electrode boundary. *Corrosion science*, 49(2):319-346, 2007.

[10] Stefan Scheiner and Christian Hellmich. Finite volume model for diffusion-and activation-controlled pitting corrosion of stainless steel. *Computer Methods in Applied Mechanics and Engineering*, 198(37-40):2898-2910, 2009.

[11] Marianne Liebi, Viviane Lutz-Bueno, Manuel Guizar-Sicairos, Bernd M. Schönbauer, Johannes Eichler, Elisabeth Martinelli, Jörg F. Löffler, Annelie Weinberg, Helga Lichteneger, and Tilman A. Grünewald. 3d nanoscale analysis of bone healing around degrading mg implants evaluated by x-ray scattering tensor tomography. *Acta Biomaterialia*, 134:804-817, October 2021.

[12] Olek C Zienkiewicz and Robert L Taylor. *The finite element method for solid and structural mechanics*. Elsevier, 2005.

[13] Antonia Torcasio, Xiaolei Zhang, Hans Van Oosterwyck, Joke Duyck, and G Harry van Lenthe. Use of micro-ct-based finite element analysis to accurately quantify peri-implant bone strains: a validation in rat tibiae. *Biomechanics and modeling in mechanobiology*, 11:743-750, 2012.

[14] Romane Blanchard, Alexander Dejaco, Evi Bongaers, and Christian Hellmich. Intravoxel bone micromechanics for microct-based finite element simulations. *Journal of biomechanics*, 46(15):2710-2721, 2013.

[15] Andrea Ojdanic, Erhard Schafler, Jelena Horky, Dmytro Orlov, and Michael Zehetbauer. Strengthening of a biodegradable mg-zn-ca alloy zx50 after processing by hpt and heat treatment. In *Magnesium Technology 2018*, pages 277-282. Springer, 2018.

[16] Heinz Parkus. *Mechanik der festen Körper [Mechanics of solids]*. Springer-Verlag, 2 edition, 1966. in German.

[17] Gillian D Muir and Ian Q Whishaw. Ground reaction forces in locomoting hemiparkinsonian rats: a definitive test for impairments and compensations. *Experimental Brain Research*, 126:307-314, 1999.

[18] Tim Wehner, Uwe Wolfram, Thomas Henzler, Frank Niemeyer, Lutz Claes, and Ulrich Simon. Internal forces and moments in the femur of the rat during gait. *Journal of Biomechanics*, 43(13):2473-2479, 2010.

[19] Ferdinand Pierre Beer, Elwood Russell Johnston, John T DeWolf, and David Francis Mazurek. * statics and Mechanics of Materials*, volume 1. mcgraw-Hill New York, 2011.

[20] MS Qatu. In-plane vibration of slightly curved laminated composite beams. *Journal of Sound and Vibration*, 159(2):327-338, 1992.

[21] Mohamad S Qatu. Theories and analyses of thin and moderately thick laminated composite curved beams. *International Journal of Solids and Structures*, 30(20):2743-2756, 1993.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

[22] Ferdinand Pierre Beer, Elwood Russell Johnston, John T DeWolf, David Francis Mazurek, and Sanjeev Sanghi. *Mechanics of materials*, volume 1. mcgraw-Hill New York, 1992.

[23] Luc Dormieux, A Molinari, and Djimedo Kondo. Micromechanical approach to the behavior of poroelastic materials. *Journal of the Mechanics and Physics of Solids*, 50(10):2203-2231, 2002.

[24] Wolfgang Kreher. Residual stresses and stored elastic energy of composites and polycrystals. *Journal of the Mechanics and Physics of Solids*, 38(1):115-128, 1990.

[25] Jean-François Barthélémy and Luc Dormieux. A micromechanical approach to the strength criterion of drucker-prager materials reinforced by rigid inclusions. *International Journal for Numerical and Analytical Methods in Geomechanics*, 28(7-8):565-582, 2004.

[26] Andreas Fritsch, Luc Dormieux, Christian Hellmich, and Julien Sanahuja. Micromechanics of crystal interfaces in polycrystalline solid phases of porous media: fundamentals and application to strength of hydroxyapatite biomaterials. *Journal of Materials Science*, 42:8824-8837, 2007.

[27] Herbert W Müllner, Andreas Fritsch, Christoph Kohlhauser, R Reihsner, C Hellmich, D Godlinski, A Rota, Robert Slesinski, and Josef Eberhardsteiner. Acoustical and poromechanical characterisation of titanium scaffolds for biomedical applications. *Strain*, 44(2):153-163, 2008.

[28] Karin Hofstetter, Ch Hellmich, Josef Eberhardsteiner, and Herbert A Mang. Micromechanical estimates for elastic limit states in wood materials, revealing nanostructural failure mechanisms. *Mechanics of Advanced Materials and Structures*, 15(6-7):474-484, 2008.

[29] Thomas K Bader, Karin Hofstetter, Christian Hellmich, and Josef Eberhardsteiner. Poromechanical scale transitions of failure stresses in wood: from the lignin to the spruce level. *ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik*, 90(10-11):750-767, 2010.

[30] Bernhard Pichler and Christian Hellmich. Upscaling quasi-brittle strength of cement paste and mortar: A multi-scale engineering mechanics model. *Cement and Concrete Research*, 41(5):467-476, 2011.

[31] Bernhard Pichler, Christian Hellmich, Josef Eberhardsteiner, Jaromír Wasserbauer, Pipat Termkhajornkit, Rémi Barbarulo, and Gilles Chanvillard. Effect of gel-space ratio and microstructure on strength of hydrating cementitious materials: An engineering micromechanics approach. *Cement and Concrete Research*, 45:55-68, 2013.

[32] Markus Königsberger, Michal Hlobil, Brice Delsaute, Stéphanie Staquet, Christian Hellmich, and Bernhard Pichler. Hydrate failure in itz governs concrete strength: A micro-to-macro validated engineering mechanics model. *Cement and Concrete Research*, 103:77-94, 2018.

[33] Jean Salençon. *Handbook of continuum mechanics: General concepts thermoelasticity*. Springer Science & Business Media, 2001.

[34] Charles H Turner. Skeletal adaptation to mechanical loading. *Clinical Reviews in Bone and Mineral Metabolism*, 5:181-194, 2007.

[35] A Michael Parfitt. The physiologic and clinical significance of bone histomorphometric data. *Bone histomorphometry: techniques and interpretation*, pages 143-223, 1983.

[36] Claire Morin, Stéphane Avril, and Christian Hellmich. Non-affine fiber kinematics in arterial mechanics: a continuum micromechanical investigation. *ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik*, 98(12):2101-2121, 2018.

[37] Claire Morin, Christian Hellmich, Zeineb Nejim, and Stéphane Avril. Fiber rearrangement and matrix compression in soft tissues: multiscale hypoelasticity and application to tendon. *Frontiers in bioengineering and biotechnology*, 9:725047, 2021.

[38] H Clarke Anderson, Rama Garimella, and Sarah E Tague. The role of matrix vesicles in growth plate development and biomineralization. *Front Biosci*, 10(1):822-837, 2005.

[39] Claire Morin and Christian Hellmich. Mineralization-driven bone tissue evolution follows
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

from fluid-to-solid phase transformations in closed thermodynamic systems. *Journal of theoretical biology*, 335:185-197, 2013.

[40] Avinash H Ambre, Dinesh R Katti, and Kalpana S Katti. Biomineralized hydroxyapatite nanoclay composite scaffolds with polycaprolactone for stem cell-based bone tissue engineering. *Journal of Biomedical Materials Research Part A*, 103(6):2077-2101, 2015.

[41] Elena Macías-Sánchez, Nadezda V Tarakina, Danail Ivanov, Stéphane Blouin, Andrea M Berzlanovich, and Peter Fratzl. Spherulitic crystal growth drives mineral deposition patterns in collagen-based materials. *Advanced Functional Materials*, 32(31):2200504, 2022.

[42] N Reznikov, B Hoac, DJ Buss, WN Addison, NMT Barros, and MD McKee. Biological stenciling of mineralization in the skeleton: local enzymatic removal of inhibitors in the extracellular matrix. *Bone*, 138:115447, 2020.

[43] K. S. Prostak and S. Lees. Visualization of crystal-matrix structure. in situ demineralization of mineralized turkey leg tendon and bone. *Calcified Tissue International*, 59(6):474-479, December 1996.

[44] Sidney Lees, Laurence C Bonar, and Herbert A Mook. A study of dense mineralized tissue by neutron diffraction. *International Journal of Biological Macromolecules*, 6(6):321-326, 1984.

[45] Christian Hellmich and F-J Ulm. Average hydroxyapatite concentration is uniform in the extracollagenous ultrastructure of mineralized tissues: evidence at the 1-10-µm scale. *Biomechanics and Modeling in Mechanobiology*, 2(1):21-36, 2003.

[46] Christian Hellmich and Franz-Josef Ulm. Micromechanical model for ultrastructural stiffness of mineralized tissues. *Journal of engineering mechanics*, 128(8):898-908, 2002.

[47] Christian Hellmich, Jean-François Barthélémy, and Luc Dormieux. Mineral-collagen interactions in elasticity of bone ultrastructure-a continuum micromechanics approach. *European Journal of Mechanics-A/Solids*, 23(5):783-810, 2004.

[48] Lukas Eberhardsteiner, Christian Hellmich, and Stefan Scheiner. Layered water in crystal interfaces as source for bone viscoelasticity: arguments from a multiscale approach. *Computer methods in biomechanics and biomedical engineering*, 17(1):48-63, 2014.

[49] Andreas Fritsch, Christian Hellmich, and Luc Dormieux. Ductile sliding between mineral crystals followed by rupture of collagen crosslinks: experimentally supported micromechanical explanation of bone strength. *Journal of theoretical biology*, 260(2):230-252, 2009.

[50] Valentina Kumbolder, Claire Morin, Stefan Scheiner, and Christian Hellmich. Hierarchical elastoplasticity of cortical bone: observations, mathematical modeling, validation. *Mechanics of Materials*, 198:105140, 2024.

[51] M Touratier. An efficient standard plate theory. *International journal of engineering science*, 29(8):901-916, 1991.

[52] Patricia Kuttke, Aleš Kurfürst, Stefan Scheiner, and Christian Hellmich. Sequential 1d/2d finite element analyses of tramway rails under bending and restrained torsion, based on the principle of virtual power. *Mechanics of Advanced Materials and Structures*, 28(11):1147-1169, 2021.

[53] Erasmo Carrera, Matteo Filippi, and Enrico Zappino. Free vibration analysis of rotating composite blades via carrera unified formulation. *Composite Structures*, 106:317-325, 2013.

[54] Alfonso Pagani and Erasmo Carrera. Large-deflection and post-buckling analyses of laminated composite beams by carrera unified formulation. *Composite Structures*, 170:40-52, 2017.

[55] Paul Germain. The method of virtual power in continuum mechanics. part 2: Microstructure. *SIAM Journal on Applied Mathematics*, 25(3):556-575, 1973.

[56] Raphael Höller, Mehdi Aminbaghai, Lukas Eberhardsteiner, Josef Eberhardsteiner, Ronald Blab, Bernhard Pichler, and Christian Hellmich. Rigorous amendment of vlasov's theory for thin elastic plates on elastic winkler foundations, based on the principle of virtual power. *European Journal of Mechanics-A/Solids*, 73:449-482, 2019.

[57] Raphael Scharf, Bernhard Pichler, Roman Heissenberger, Bernd Moritz, and Christian
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

Hellmich. Data-driven analytical mechanics of aging viscoelastic shotcrete tunnel shells.
*Acta Mechanica*, 233(8):2989-3019, 2022.

[58] Danilo Capecchi. History of virtual work laws: a history of mechanics prospective, volume 42. Springer Science & Business Media, 2012.

[59] Francesco dell'Isola and Anil Misra. Principle of virtual work as foundational framework for metamaterial discovery and rational design. *Comptes Rendus. Mécanique*, 351(S3):1-25, 2023.