17 Warum Knochen sich umbauen: FEM mit Rückkopplung

Kapitel 16 gab jedem Dreieck eine feste Dichte und daraus einen festen E-Modul. Die FEM antwortete mit Verschiebungen und Spannungen, aber das Material hörte nicht zu. Jetzt schließen wir den Kreis:

\[ \rho \longrightarrow E(\rho) \longrightarrow \text{FEM} \longrightarrow \text{mechanischer Reiz} \longrightarrow \rho_{\mathrm{neu}}. \]

Das ist eine neue Art von Rechnung. Nicht der Löser ist neu, sondern die Rückkopplung: Das Ergebnis eines statischen FEM-Laufs verändert die Materialkarte für den nächsten Lauf. Wir untersuchen damit ein rechteckiges, spongiöses Ersatzfenster. Es entstehen dichtere und weniger dichte Bereiche, aber weder Löcher noch einzelne Knochenbälkchen.

Naheliegende Vermutung

Vermutung: „Wo die Spannung groß ist, wächst Knochen. Man wiederholt das einfach, bis ein Trabekelmuster entsteht.“

Was stattdessen stimmt: Eine rein lokale Regel kann benachbarte Elemente abwechselnd steif und weich machen. Das Schachbrett ist ein Netzartefakt, kein Trabekelmuster. Außerdem enthält unser Modell weder Zellen noch Stoffwechsel oder reale Umbauzeiten. Es ist ein prüfbares mechanisches Ersatzmodell für Dichteanpassung.

Lernziele

Nach diesem Kapitel kannst du …

… die Remodeling-Schleife als Folge statischer FEM-Rechnungen erklären,
… aus Spannung und Dehnung die Formänderungsenergiedichte und daraus einen Reiz je Knochenmasse berechnen,
… Abbau, Totzone und Aufbau in einer begrenzten Dichte-Regel unterscheiden,
… erklären, warum ein lokales Update Schachbrettmuster erzeugen kann und wie eine nichtlokale Mittelung dagegen wirkt,
… Stationarität, Homöostase und fehlende Konvergenz auseinanderhalten,
… die Grenzen eines isotropen Dichtekontinuums ehrlich benennen.

17.1 Viele statische Rechnungen — keine biologische Uhr

Das Ersatzfenster ist 60 mm lang, 30 mm hoch und 10 mm tief. Links halten wir alle Knoten fest. Rechts wirken zwei getrennte, kraftgesteuerte Lastfälle:

Lastfall	Resultierende	Gewicht
Biegung	600 N nach unten	0,65
Zug	600 N nach rechts	0,35

Die Gewichte addieren sich zu eins. Sie sind eine offengelegte Mischung zweier didaktischer Lastpfade, keine gemessenen Häufigkeiten von Gehen oder Treppensteigen. In jeder Recheniteration lösen wir beide Lastfälle mit denselben Kräften und Lagern. Danach erst ändern wir die Dichten synchron.

Iteration ist nicht Zeit

Iteration 40 bedeutet weder Tag 40 noch Monat 40. Die Schrittweite steuert eine numerische Fixpunktiteration. Ohne kalibriertes biologisches Zeitgesetz darf ihre Nummer nicht in Lebenszeit übersetzt werden.

Code

# Der Modellvertrag in einer kleinen, sichtbaren Tabelle.
parameter = [
    ("Geometrie", "60 x 30 x 10 mm"),
    ("Lehrnetz", "45 Knoten, 64 P1-Dreiecke"),
    ("Dichtebereich", "0,2 bis 1,0 g/cm3"),
    ("Startdichte", "0,5 g/cm3"),
    ("Lastgewichte", "0,65 / 0,35"),
]
for name, wert in parameter:
    print("%-16s %s" % (name, wert))

Geometrie        60 x 30 x 10 mm
Lehrnetz         45 Knoten, 64 P1-Dreiecke
Dichtebereich    0,2 bis 1,0 g/cm3
Startdichte      0,5 g/cm3
Lastgewichte     0,65 / 0,35

17.2 Ein Dreieck speichert Verformungsenergie

Drückst du eine Feder langsam zusammen, steckt anschließend Energie in ihrer Verformung. Ein Dreieck speichert ebenso elastische Energie. Pro Volumen heißt sie Formänderungsenergiedichte \(U\). Mit der Engineering-Voigt-Schreibweise

\[ \boldsymbol\varepsilon= [\varepsilon_{xx},\varepsilon_{yy},\gamma_{xy}]^T, \qquad \boldsymbol\sigma=[\sigma_{xx},\sigma_{yy},\tau_{xy}]^T \]

gilt

\[ U=\frac12\,\boldsymbol\sigma^T\boldsymbol\varepsilon. \]

Der Schubterm ist damit \(\tau_{xy}\gamma_{xy}\). Das ist wichtig: \(\gamma_{xy}\) ist die Engineering-Schubdehnung, nicht die tensorielle Komponente \(\varepsilon_{xy}=\gamma_{xy}/2\).

Nehmen wir für ein Dreieck \(\boldsymbol\sigma=[2{,}0;\,1{,}0;\,0{,}5]\) Pa und \(\boldsymbol\varepsilon=[2{,}0;\,1{,}0;\,1{,}0]\). Dann ist

\[ U=\tfrac12(2\cdot2+1\cdot1+0{,}5\cdot1)=2{,}75\ \mathrm{J/m^3}. \]

Code

import os, sys
sys.path.insert(0, os.path.join("..", "programme", "gemeinsam"))
from remodeling import formenergie_dreieck

spannung = [2.0, 1.0, 0.5]
dehnung = [2.0, 1.0, 1.0]
u = formenergie_dreieck(spannung, dehnung)
print("U =", u, "J/m3")

U = 2.75 J/m3

17.3 Energie je Masse wird zum Reiz

Dieselbe Energie pro Volumen bedeutet bei wenig vorhandener Knochenmasse mehr Energie je Masse. Für Dreieck \(i\) und Lastfall \(l\) bilden wir deshalb

\[ S_i=\frac{\sum_l p_l U_{i,l}}{\rho_{\mathrm m,i}}, \qquad \rho_{\mathrm m,i}=1000\,\rho_i. \]

Im Morgan-Gesetz steht \(\rho_i\) in g/cm³; nach Multiplikation mit 1000 steht \(\rho_{\mathrm m,i}\) in kg/m³. Weil \(U\) die Einheit J/m³ besitzt, hat \(S\) die Einheit J/kg. Energie je Knochenmasse wird in adaptiven Knochenmodellen als mechanische Rückkopplungsgröße verwendet (Huiskes u. a. 1987; Weinans u. a. 1992). Unsere lineare Lastfallmischung ist eine eigene didaktische Modellentscheidung.

Beispiel: \(U=1{,}0\cdot10^4\) J/m³ und \(\rho=0{,}5\) g/cm³ ergeben \(\rho_\mathrm m=500\) kg/m³ und \(S=20\) J/kg.

Code

energie = 1.0e4       # J/m3
dichte = 0.5          # g/cm3
massendichte = 1000.0 * dichte
reiz = energie / massendichte
print("Massendichte =", massendichte, "kg/m3")
print("Reiz         =", reiz, "J/kg")

Massendichte = 500.0 kg/m3
Reiz         = 20.0 J/kg

Für den Zielreiz brauchen wir einen festen Maßstab. Wir rechnen einmal einen gleichmäßigen Referenzzustand mit \(\rho_\mathrm{ref}=0{,}5\) g/cm³ auf einem 16×8-Netz und mitteln sein Reizfeld flächengewichtet:

\[ S_0=\frac{\sum_i A_iS_i^\mathrm{ref}}{\sum_i A_i} =19{,}4294626265\ \mathrm{J/kg}. \]

\(S_0\) bleibt danach fest — auch bei anderer Startdichte, Schrittweite, Netzfeinheit oder Lastmischung. Er ist ein Referenzwert dieses Modells, kein universell gemessener Sollwert für Knochen.

17.4 Der schlechte erste Versuch: jedes Element hört nur sich selbst

Die verlockende Regel lautet: Ist \(S_i>S_0\), erhöhe die Dichte; ist \(S_i<S_0\), verringere sie. Unter Kraftlast wird ein etwas weicheres Element stärker verformt. Sein Reiz steigt, sein Update kann es wieder versteifen — und zugleich Last in ein Nachbarelement umlenken. Auf einem diskreten Netz kann diese Rückkopplung benachbarte Dreiecke gegeneinander treiben.

Das Ergebnis ist oft ein Schachbrettmuster: große Sprünge auf der Skala eines einzelnen Elements. Es kann stationär werden und trotzdem falsch sein. Wie in Kapitel 12 gilt: Ein buntes Bild und ein kleiner Abbruchfehler sind noch keine Netzunabhängigkeit.

17.5 Die Nachbarschaft hört mit

Wir ersetzen den lokalen Reiz durch einen Mittelwert über die Schwerpunkte der Nachbarelemente. Ein nahes Dreieck zählt stark, ein Dreieck am Filterrand gar nicht. Das Gewicht ist ein Zeltdach über dem Abstand:

\[ \bar S_i= \begin{cases} S_i,&R=0,\\[1mm] \displaystyle\frac{\sum_j w_{ij}A_jS_j}{\sum_jw_{ij}A_j},&R>0, \end{cases} \qquad w_{ij}=\max\!\left(0,1-\frac{d_{ij}}R\right). \]

Die Fläche \(A_j\) verhindert, dass viele kleine Dreiecke mehr zählen als wenige große. \(R=0\) ist ausdrücklich der ungefilterte Kontrollfall \(\bar S_i=S_i\); ein negativer Radius ist ungültig. Für \(R>0\) trägt das Zielelement selbst immer mit Gewicht eins bei.

Der verwendete Schwerpunkt-/Flächenfilter ist eine bewusst einfache Adaption nichtlokaler Mittelung gegen Konvergenz- und Schachbrettprobleme (Calvo-Gallego u. a. 2021), nicht die identische Reproduktion des dortigen Algorithmus. Im Lehrmodell gilt \(R=30\) mm als feste physische Länge. Auf dem 8×4-Netz tragen je nach Randlage 28 bis 59 Dreiecke bei; der Median ist 43,5.

Code

from remodeling import baue_filtergewichte, glaette_reize

schwerpunkte = [(0.0, 0.0), (1.0, 0.0), (2.0, 0.0)]
flaechen = [1.0, 2.0, 1.0]
lokaler_reiz = [10.0, 20.0, 40.0]

filter_aus = baue_filtergewichte(schwerpunkte, flaechen, 0.0)
filter_an = baue_filtergewichte(schwerpunkte, flaechen, 1.5)
print("R = 0  :", glaette_reize(lokaler_reiz, filter_aus))
print("R = 1,5:", glaette_reize(lokaler_reiz, filter_an))

R = 0  : [10.0, 20.0, 40.0]
R = 1,5: [14.0, 21.25, 32.0]

Code

import matplotlib.pyplot as plt
import numpy as np

fig, achsen = plt.subplots(1, 2, figsize=(10.0, 3.8))
achsen[0].axis("off")
punkte = [(0.10, 0.58, "Dichte ρ"), (0.30, 0.82, "E(ρ)"),
          (0.60, 0.82, "FEM"), (0.80, 0.58, "Reiz S"),
          (0.46, 0.28, "Update")]
for x, y, text in punkte:
    achsen[0].text(x, y, text, ha="center", va="center", fontsize=11,
                   bbox=dict(boxstyle="round,pad=0.35", fc="#eaf2f8", ec="#2874a6"))
for start, ende in [(punkte[0], punkte[1]), (punkte[1], punkte[2]),
                    (punkte[2], punkte[3]), (punkte[3], punkte[4]),
                    (punkte[4], punkte[0])]:
    achsen[0].annotate("", xy=(ende[0], ende[1]), xytext=(start[0], start[1]),
                      arrowprops=dict(arrowstyle="->", lw=1.5, color="#34495e"))
achsen[0].set_title("Rückkopplung statt Einweg")

a = np.linspace(-0.55, 0.55, 400)
d = 0.08
g = np.where(a < -d, a + d, np.where(a > d, a - d, 0.0))
achsen[1].plot(a, g, lw=2.2)
achsen[1].axvspan(-d, d, color="tab:green", alpha=0.15, label="Totzone")
achsen[1].axhline(0.0, color="0.3", lw=0.8)
achsen[1].axvline(0.0, color="0.3", lw=0.8)
achsen[1].text(-0.36, -0.20, "Abbau", color="tab:blue")
achsen[1].text(0.26, 0.20, "Aufbau", color="tab:red")
achsen[1].set_xlabel("Antrieb a = S̄/S₀ − 1")
achsen[1].set_ylabel("Antwort g(a)")
achsen[1].legend()
achsen[1].grid(alpha=0.2)
plt.tight_layout()
plt.show()

Abbildung 17.1: Die geschlossene Remodeling-Schleife (links) und die Antwortfunktion (rechts). Innerhalb der Totzone bleibt die Dichte unverändert; außerhalb reagiert sie stetig mit Abbau oder Aufbau.

17.6 Aufbau, Abbau und die Totzone

Frosts Mechanostat liefert einen biologischen Rahmen dafür, dass Knochen auf unterschiedliche mechanische Beanspruchung unterschiedlich reagiert (Frost 2003). Wir übernehmen daraus keine scheinpräzisen biologischen Schwellen. Stattdessen definieren wir den dimensionslosen Antrieb

\[ a_i^n=\frac{\bar S_i^n}{S_0}-1 \]

und eine transparente numerische Antwort mit Totzonenbreite \(\delta=0{,}08\):

\[ g(a)= \begin{cases} a+\delta,&a<-\delta,\\ 0,&|a|\leq\delta,\\ a-\delta,&a>\delta. \end{cases} \]

Mit \(\eta=0{,}08\) g/cm³ je Recheniteration und \(\Delta\rho_{\max}=0{,}04\) g/cm³ je Recheniteration folgt

\[ \Delta\rho_i^n=\operatorname{clip} \left(\eta g(a_i^n),-\Delta\rho_{\max},\Delta\rho_{\max}\right), \]

\[ \rho_i^{n+1}=\operatorname{clip} \left(\rho_i^n+\Delta\rho_i^n,0{,}2,1{,}0\right). \]

clip bedeutet nur: Werte unter der Untergrenze werden auf die Untergrenze, Werte über der Obergrenze auf die Obergrenze gesetzt. Alle neuen Dichten werden aus derselben unveränderten Kopie von \(\rho^n\) berechnet. Ein Update mitten in der Elementschleife wäre ein anderes, von der Nummerierung abhängiges Verfahren.

Code

from remodeling import aktualisiere_dichten

dichten = [0.50, 0.50, 0.50]
reize = [0.80 * 19.4294626265,
         1.04 * 19.4294626265,
         1.30 * 19.4294626265]
neu, antriebe, aenderungen = aktualisiere_dichten(
    dichten, reize, 19.4294626265, 0.08, 0.08, 0.04, 0.2, 1.0)
for i in range(3):
    print("a=%5.2f  delta_rho=%6.3f  rho_neu=%5.3f"
          % (antriebe[i], aenderungen[i], neu[i]))

a=-0.20  delta_rho=-0.010  rho_neu=0.490
a= 0.04  delta_rho= 0.000  rho_neu=0.500
a= 0.30  delta_rho= 0.018  rho_neu=0.518

17.7 Die vollständige Schleife

Eine Iteration besteht nun aus sechs klaren Schritten:

Aus jeder Dichte mit dem Morgan-Gesetz den E-Modul bilden.
Beide festen Lastfälle lösen.
Je Dreieck \(U\) und den gewichteten lokalen Reiz \(S\) berechnen.
\(S\) über die feste physische Nachbarschaft zu \(\bar S\) mitteln.
Alle Dichten synchron und begrenzt aktualisieren.
Änderung, Zwei-Schritt-Abstand, Gleichgewicht und Diagnosen prüfen.

Das eigenständige Kapitelprogramm führt genau diese Schleife aus:

Code

import os, sys
sys.path.insert(0, os.path.join("..", "programme", "kap17"))
import knochen_remodeling as kr

ergebnis = kr.simuliere()
lokal = kr.simuliere(filter_radius=0.0)
diag = ergebnis["diagnosen"]
print("Status             :", ergebnis["status"])
print("Iterationen        :", ergebnis["iterationen"])
print("Dichtemittel       : %.5f g/cm3" % diag["dichtemittel"])
print("Masse              : %.6f kg" % diag["masse"])
print("Nachgiebigkeit     : %.6e J" % diag["nachgiebigkeit"])
print("Grenz-/Aussenanteil: %.1f %% / %.1f %%"
      % (100.0 * diag["grenzanteil"], 100.0 * diag["aussenanteil"]))
print("q_hf lokal/gefiltert: %.5f / %.5f" % (lokal["q_hf"], ergebnis["q_hf"]))

Status             : stationaer
Iterationen        : 88
Dichtemittel       : 0.42971 g/cm3
Masse              : 0.007735 kg
Nachgiebigkeit     : 3.164353e-01 J
Grenz-/Aussenanteil: 0.0 % / 3.1 %
q_hf lokal/gefiltert: 0.36135 / 0.02073

Code

import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation

xs = []
ys = []
tris = []
for punkt in ergebnis["knoten"]:
    xs.append(1000.0 * punkt[0])
    ys.append(1000.0 * punkt[1])
for dreieck in ergebnis["dreiecke"]:
    tris.append(list(dreieck))
tri = Triangulation(xs, ys, tris)

fig, achsen = plt.subplots(2, 2, figsize=(10.0, 6.6))
p0 = achsen[0, 0].tripcolor(tri, facecolors=lokal["dichten"], cmap="viridis",
                            vmin=0.2, vmax=1.0, edgecolors="0.75", linewidth=0.25)
achsen[0, 0].set_title("lokal: R = 0")
fig.colorbar(p0, ax=achsen[0, 0], label="Dichte (g/cm³)")
p1 = achsen[0, 1].tripcolor(tri, facecolors=ergebnis["dichten"], cmap="viridis",
                            vmin=0.2, vmax=1.0, edgecolors="0.75", linewidth=0.25)
achsen[0, 1].set_title("nichtlokal: R = 30 mm")
fig.colorbar(p1, ax=achsen[0, 1], label="Dichte (g/cm³)")

it = []
mittel = []
for i in range(len(ergebnis["diagnoseverlauf"])):
    it.append(i + 1)
    mittel.append(ergebnis["diagnoseverlauf"][i]["dichtemittel"])
achsen[1, 0].plot(it, mittel, color="tab:blue")
achsen[1, 0].set_ylabel("mittlere Dichte (g/cm³)")
achsen[1, 0].set_xlabel("Recheniteration")
achsen[1, 0].grid(alpha=0.2)
achsen[1, 1].semilogy(it, ergebnis["aenderungsverlauf"], label="r∞")
achsen[1, 1].semilogy(it[1:], ergebnis["zweischrittverlauf"][1:], label="r₂")
achsen[1, 1].axhline(8.0e-5, color="tab:red", ls="--", label="Toleranz")
achsen[1, 1].set_xlabel("Recheniteration")
achsen[1, 1].set_ylabel("Dichteänderung (g/cm³)")
achsen[1, 1].legend()
achsen[1, 1].grid(alpha=0.2)
for ax in achsen[0]:
    ax.set_aspect("equal")
    ax.set_xlabel("x (mm)")
    ax.set_ylabel("y (mm)")
plt.tight_layout()
plt.show()

Abbildung 17.2: Das lokale Modell (oben links) wird zwar stationär, zeigt aber ein elementgroßes Schachbrett. Der feste nichtlokale Filter (oben rechts) erhält den großskaligen Lastpfad und senkt den Hochfrequenzindikator um 94,3 %. Unten: Dichtemittel sowie Ein- und Zwei-Schritt-Änderung des gefilterten Laufs.

Das lokale Feld wirkt wie ein Dreiecksmosaik. Der Hochfrequenzindikator

\[ q_\mathrm{hf}=\sqrt{\frac{\sum_i A_i (\rho_i-\langle\rho\rangle_{\mathcal N_i})^2}{\sum_iA_i}} \]

misst die Abweichung vom Mittel der Kantennachbarn. Er fällt von 0,36135 auf 0,02073 g/cm³ — um 94,3 %. Ein kleinerer Kantensprung allein wäre kein Beweis, denn auch ein sinnvoller glatter Gradient besitzt Kantensprünge.

Code

from IPython.display import HTML
from matplotlib.animation import FuncAnimation
import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation

fig, ax = plt.subplots(figsize=(7.5, 3.6))
tri_anim = Triangulation(xs, ys, tris)
bild = ax.tripcolor(tri_anim, facecolors=ergebnis["dichteverlauf"][0],
                    cmap="viridis", vmin=0.2, vmax=1.0,
                    edgecolors="0.75", linewidth=0.25)
fig.colorbar(bild, ax=ax, label="Dichte (g/cm³)")
ax.set_aspect("equal")
ax.set_xlabel("x (mm)")
ax.set_ylabel("y (mm)")

def zeichne_remodeling(nummer):
    bild.set_array(ergebnis["dichteverlauf"][nummer])
    ax.set_title("Recheniteration %d — keine Zeitangabe" % nummer)
    return bild,

bilder = []
for nummer in range(0, ergebnis["iterationen"] + 1, 3):
    bilder.append(nummer)
if bilder[-1] != ergebnis["iterationen"]:
    bilder.append(ergebnis["iterationen"])
ani = FuncAnimation(fig, zeichne_remodeling, frames=bilder,
                    interval=180, blit=False)
plt.close(fig)
HTML(ani.to_jshtml())

Im HTML lässt sich die Folge vor- und zurückspulen. Die feste Farbskala macht sichtbar, dass zuerst die Materialmenge abnimmt und sich danach der glatte räumliche Verlauf einpendelt. Für den Druck zeigt die nächste Abbildung sechs beschriftete Zustände derselben Rechnung.

Code

import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation

auswahl = [0, 5, 10, 20, 40, ergebnis["iterationen"]]
fig, achsen = plt.subplots(2, 3, figsize=(9.5, 5.5), sharex=True, sharey=True)
achsen_liste = []
for zeile in achsen:
    for ax in zeile:
        achsen_liste.append(ax)
for ax, nummer in zip(achsen_liste, auswahl):
    p = ax.tripcolor(tri, facecolors=ergebnis["dichteverlauf"][nummer],
                     cmap="viridis", vmin=0.2, vmax=1.0,
                     edgecolors="0.75", linewidth=0.2)
    ax.set_title("Iteration %d" % nummer)
    ax.set_aspect("equal")
    ax.set_xlabel("x (mm)")
achsen[0, 0].set_ylabel("y (mm)")
achsen[1, 0].set_ylabel("y (mm)")
fig.colorbar(p, ax=achsen_liste, shrink=0.78, label="Dichte (g/cm³)")
plt.show()

Abbildung 17.3: Sechs Zustände derselben Recheniteration mit fester Farbskala. Die Startkarte ist gleichmäßig; danach bildet sich ein glatter großskaliger Lastpfad. Die Nummer ist eine Recheniteration, keine Zeitangabe.

17.8 Wann ist die Schleife fertig?

Wir messen nach jedem Update

\[ r_\infty^n=\max_i|\rho_i^{n+1}-\rho_i^n|. \]

Der Lauf heißt stationär, wenn \(r_\infty<8\cdot10^{-5}\) g/cm³ in drei aufeinanderfolgenden Iterationen gilt. Nach spätestens 200 Iterationen endet er sonst als „nicht konvergiert“. Zusätzlich beobachten wir \(r_2^n=\max_i|\rho_i^{n+1}-\rho_i^{n-1}|\): Wird nur dieser Wert klein, pendelt das Verfahren zwischen zwei Feldern und heißt „oszillierend“.

Stationär bedeutet aber nicht automatisch homöostatisch. Ein Element kann an einer Dichtegrenze festhängen, obwohl sein Reiz außerhalb der Totzone liegt. Deshalb geben wir getrennt den Flächenanteil an Dichtegrenzen und den Anteil außerhalb der Totzone aus. Nur stationär und ohne äußeren Anteil nennen wir in diesem didaktischen Modell homöostatisch. Der Referenzlauf ist stationär nach 88 Iterationen, aber mit 3,1 % außerhalb der Totzone nicht vollständig homöostatisch.

Weitere Kontrollgrößen sind die Masse \(m=t\sum_iA_i\rho_{\mathrm m,i}\), die gewichtete Nachgiebigkeit \(C=\sum_lp_l\mathbf f_l^T\mathbf u_l\) sowie Kraft- und Momentengleichgewicht. Beim Referenzlauf beträgt die Masse 0,007735 kg und \(C=0{,}3164353\) J.

17.9 Gegen eine unabhängige Rechnung geprüft

Der reine Python-Code wurde nicht nur gegen sich selbst getestet. Eine unabhängige NumPy/scikit-fem-Implementierung rechnet dieselben Netze und Lasten. Die Freigabe umfasst harte Assertions; ein Fehler beendet das Validierungsprogramm mit einem Fehlercode.

Prüfung	Ergebnis
Verschiebung, Spannung, Dehnung, Reaktion	relative Fehler < \(1{,}1\cdot10^{-13}\)
Reizfeld	relativer Fehler \(2{,}0\cdot10^{-13}\)
vollständiger synchroner Update-Schritt	relativer Fehler \(2{,}7\cdot10^{-14}\)
Kraft-/Momentengleichgewicht	relative Fehler < \(4{,}5\cdot10^{-14}\)
12×6 gegen 16×8, Enddichtefeld	3,28 %
zweite Diagonalenorientierung	3,12 %
halbe Schrittweite, Enddichtefeld	0,51 %
Startdichten 0,4 / 0,55 g/cm³	3,43 % / 3,87 %
\(q_\mathrm{hf}\), lokal → nichtlokal	0,36135 → 0,02073

Alle festgelegten Grenzen werden eingehalten. Die Details und das maschinenlesbare Ergebnis liegen in validierung/kap17_remodeling/. Besonders wichtig: Drei Netze werden stationär, und zwischen den zwei feinsten Netzen bleibt nicht nur ein Mittelwert, sondern das auf ein gemeinsames Raster projizierte Feld innerhalb von 5 %.

17.10 Probier es aus: die Remodeling-Werkstatt

Die folgende Browserrechnung benutzt dasselbe 45-Knoten-Netz, dieselben zwei Lastfälle und dasselbe feste \(S_0\). Verändere Filterradius, Totzone oder das Gewicht des Biegelastfalls. Es werden bewusst nur 15 Schritte gezeigt. Beim ersten Start muss der Browser Pyodide und Matplotlib laden; auf der Referenzhardware dauerte dieser Kaltstart 55,7 s. Danach benötigten fünf Wiederholungen im Median 1,41 s (maximal 1,56 s).

Vorhersage-Punkt

Setze zuerst FILTER_RADIUS = 0.0. Erwartest du ein glatteres oder ein stärker wechselndes Dichtefeld? Begründe deine Vorhersage, bevor du ausführst.

Abbildung 17.4: Vorgerenderte Ausgangseinstellung der Remodeling-Werkstatt: 15 Recheniterationen mit R = 30 mm, Totzone 0,08 und Biegeanteil 0,65. Im HTML ersetzt die Browserrechnung dieses Bild durch dein Ergebnis.

17.11 Übungen

17.11.1 Übung 17.1 — Verstehen

Für \(delta=0{,}08\) entscheide ohne Programm, ob abgebaut, nichts geändert oder aufgebaut wird:

\(\bar S/S_0=0{,}75\),
\(\bar S/S_0=1{,}05\),
\(\bar S/S_0=1{,}30\).

Berechne jeweils \(a\) und das Vorzeichen von \(g(a)\).

Lösung zu Übung 17.1

\(a=-0{,}25< -0{,}08\): Abbau, \(g=-0{,}17\).
\(a=0{,}05\) liegt in der Totzone: keine Änderung, \(g=0\).
\(a=0{,}30>0{,}08\): Aufbau, \(g=0{,}22\).

17.11.2 Übung 17.2 — Verändern

Halbiere und verdopple \(eta\). Vergleiche Status, Iterationszahl und Endfeld. Setze danach \(R=0\) und vergleiche den Hochfrequenzindikator. Warum ist ein schneller stationärer Lauf nicht automatisch der bessere Lauf?

Lösung zu Übung 17.2

Das Skript loesungen/kap17_ue2.py rechnet alle vier Fälle. Mit \(\eta=0{,}04/0{,}08/0{,}16\) werden sie nach 149/88/69 Iterationen stationär. Der lokale Fall ist bereits nach 48 Iterationen stationär, hat aber \(q_\mathrm{hf}=0{,}36135\) statt 0,02073. Er erfüllt das Abbruchkriterium, nicht die Artefaktprüfung.

17.11.3 Übung 17.3 — Übertragen

Ergänze als dritten Lastfall 600 N diagonal nach rechts unten. Gib ihm Gewicht 0,20 und skaliere die bisherigen Gewichte auf 0,52 und 0,28. Sage vor dem Lauf voraus, wie sich der Lastpfad ändert. Interpretiere das Ergebnis nur für diese gewählte Lastmischung.

Lösung zu Übung 17.3

loesungen/kap17_ue3.py nutzt die öffentliche lastfaelle-Schnittstelle. Der erweiterte Lauf wird nach 132 Iterationen stationär; sein Dichtemittel beträgt 0,42273 g/cm³. Die diagonale Kraft koppelt Zug und Biegung im selben Lastfall. Das Resultat ist keine medizinische Vorhersage.

17.12 Was das Modell zeigt — und was nicht

Das Modell zeigt nachvollziehbar, wie Rückkopplung, Totzone, Begrenzung und räumliche Mittelung aus wiederholten FEM-Lösungen eine Dichtekarte formen. Es zeigt ebenso, dass Konvergenz mehrere Prüfungen braucht: Updatefehler, Zwei-Schritt-Abstand, Gleichgewicht, Grenzanteile, Netz- und Parameterstudien.

Es zeigt keine echten Trabekel. Jedes Dreieck bleibt ein isotropes, linear-elastisches Dichtekontinuum; Poren werden nur über die Scheindichte vertreten. Das Modell enthält keine Osteoblasten, Osteoklasten, Hormone, Ermüdung, Schädigung, Fraktur oder Heilung. Der euklidische Filter ist nur für das konvexe Rechteck geprüft und dürfte bei Löchern oder Spalten nicht ungeprüft über getrennte Flächen hinweg mitteln. Die Lastfälle sind didaktisch, die Geometrie nicht patientenspezifisch und die Iteration keine biologische Zeit.

Die wichtigste neue Einsicht ist deshalb nicht ein bestimmtes Farbmuster, sondern die geschlossene Rechenkette: Ein Material kann auf die mechanische Antwort seines eigenen FEM-Modells reagieren — vorausgesetzt, wir formulieren und prüfen diese Rückkopplung ebenso streng wie die FEM selbst.

Calvo-Gallego, José Luis, Peter Pivonka, José Manuel García-Aznar, und Javier Martínez-Reina. 2021. „A novel algorithm to resolve lack of convergence and checkerboard instability in bone adaptation simulations using non-local averaging“. International Journal for Numerical Methods in Biomedical Engineering 37 (2): e3419. https://doi.org/10.1002/cnm.3419.

Frost, Harold M. 2003. „Bone’s mechanostat: a 2003 update“. The Anatomical Record Part A 275A (2): 1081–101. https://doi.org/10.1002/ar.a.10119.

Huiskes, Rik, H. Weinans, H. J. Grootenboer, M. Dalstra, B. Fudala, und T. J. Slooff. 1987. „Adaptive bone-remodeling theory applied to prosthetic-design analysis“. Journal of Biomechanics 20 (11–12): 1135–50. https://doi.org/10.1016/0021-9290(87)90030-3.

Weinans, H., R. Huiskes, und H. J. Grootenboer. 1992. „The behavior of adaptive bone-remodeling simulation models“. Journal of Biomechanics 25 (12): 1425–41. https://doi.org/10.1016/0021-9290(92)90056-7.