Am Ende von 1 stand eine Einladung: Nimm fem_balken.py, das Programm, das den brennenden Träger von der Wärme bis zur Verformung durchrechnet, und mach es zu deinem — ein anderes Bauteil, ein anderes Material, eine andere Frage. Dieses Kapitel nimmt die Einladung beim Wort. Wir wechseln das Reich: vom Stahlträger unter der Betondecke zu einem Oberschenkelknochen unter dem Gewicht eines Menschen. Und wir stellen fest, dass die Maschine unverändert bleibt.
Das ist der eigentliche Grund, dieses Kapitel zu schreiben. Fünfzehn Kapitel lang sah es aus, als ginge es um einen Träger im Feuer. In Wahrheit ging es die ganze Zeit um eine Methode. Der Beweis dafür ist, dass wir sie ohne eine einzige neue Löser-Zeile auf etwas Lebendiges übertragen können. Ein Knochen ist kein Stahl — aber für kleine, kurze Belastungen ist er ein Träger wie jeder andere: Er hat ein Netz, eine Last, ein Auflager und Spannungen, die wir lesen können.
Beim Gehen erreicht die Kraft, mit der der Hüftkopf in die Pfanne drückt, in Messungen mit instrumentierten Prothesen grob das Zweieinhalbfache des Körpergewichts (Bergmann u. a. 2001) — bei einem Menschen von 75 kg also rund 1,8 kN. Die feinen Knochenbälkchen im Inneren des Knochens sind nicht zufällig angeordnet; Belastung ist einer von mehreren Einflüssen auf ihre Architektur. Dein Programm aus 1 kann einen bewusst vereinfachten Ausschnitt davon untersuchen.
HinweisEin neues Reich — kein neuer Erzählstrang
Die Brand-Szene ruht ab hier. Dieses Kapitel ist eine Zugabe nach dem Finale: Wir feiern, dass die Methode wandert, statt den Träger weiter im Feuer zu lassen. Der Wärmeteil aus 1 bleibt stehen, wo er ist — der Knochen wird rein mechanisch gerechnet.
Lernziele
Nach diesem Kapitel kannst du …
… erklären, warum ein tragender Knochen dasselbe FEM-Problem ist wie ein Träger — Netz, Last, Auflager, Spannungen — und warum die Hüftkontaktkraft am Gelenk wirkt, der Muskel aber nicht als zweite Last zusätzlich angesetzt wird,
… begründen, warum das Material eines Knochens ortsabhängig ist, und die Kette kalibrierter CT-Wert → Scheindichte → E-Modul erklären, ohne ein Röntgenbild, einen CT-Wert und eine Dichte gleichzusetzen,
… die Assemblierung mit elementweise verschiedenem E ausführen — den Baustein aus Kapitel 14, jetzt mit einer Materialtabelle je Dreieck,
… E-Modul-, Dehnungs- und Spannungsfelder unterscheiden und das Ergebnis mit der Prüf-Checkliste aus Kapitel 12 kontrollieren,
… das Wolffsche Gesetz und moderne Remodeling-Modelle als Ausblick auf eine Rückkopplung erklären, ohne daraus schon einen Algorithmus abzuleiten,
… die bewussten Grenzen benennen (2D statt 3D, isotrop statt faserig, statische Ersatzrechnung statt Wachstum) — Ehrlichkeit wie in Kapitel 12 und
WarnungNaheliegende Vermutung
Vermutung:„Biologie ist zu weich und lebendig für die harte Mathematik der FEM.”
Warum sie naheliegt: Knochen wächst, heilt, verändert sich — nichts davon kennt der starre Stahlträger. Und „weich” klingt nach etwas, das man nicht wie einen Balken durchrechnen kann.
Was stattdessen stimmt: Für kleine, kurzzeitige Verformungen lässt sich ein Knochen näherungsweise linear-elastisch rechnen — mit denselben Gleichungen wie der Träger. Die harte äußere Schale (Kompakta) liegt mit etwa 17 GPa zwar weit unter Stahl (210 GPa), ist aber keineswegs weich wie Gewebe; sie ist steifer als viele Kunststoffe. Das Modell ist bewusst linear und isotrop. Dass sich Knochen langfristig an Belastung anpassen, öffnet am Kapitelende eine Tür — es ersetzt aber nicht die statische Rechnung, sondern baut auf ihr auf.
16.1 Der Transfer, ausdrücklich
Legen wir die beiden Bauteile nebeneinander. Der Träger aus Kapitel 14 war ein Rechteck aus Stahl, unten vom Brand erwärmt, an einem Ende eingespannt. Der Femurschnitt ist eine krumme Kontur aus Knochen, oben von der Hüftkraft belastet, am Schaftfuß festgehalten. Verschieden sieht das aus — gleich ist die Frage: Wie verformt sich das Bauteil, und wo sitzen die Spannungen?
Zwei Dinge übersetzen wir, eine Zutat ist wirklich neu.
Übersetzung 1 — die Last kommt vom Gelenk. Das Hüftgelenk überträgt die Last auf den Hüftkopf. Die gemessene Hüftkontaktkraft von 1,8 kN ist bereits die Resultierende aus Körpergewicht und Muskelzug — sie enthält die Muskelwirkung schon. Wir setzen deshalb keine zweite, erfundene Muskelkraft zusätzlich an; das würde die Wirkung doppelt zählen. Das Auflager ist der abgeschnittene Schaftfuß: eine Modell-Lagerung, kein echtes Gelenk. Ohne Wärmeteil entfällt die thermische rechte Seite von Kapitel 14 vollständig.
Übersetzung 2 — Misstrauen als Handwerk. Die Checkliste aus Kapitel 12 („Traue keinem bunten Bild”) wird zur biologischen Frage: Ist eine Spannungsspitze am Schenkelhals echt, oder ist sie nur Kerbwirkung an einer zu grob modellierten Ecke? Wir werden genau das prüfen.
Die eine neue Zutat — das Material ist eine Landkarte. Der Knochen hat kein einheitliches Material. Sein E-Modul ist ortsabhängig: außen die steife kortikale Schale, innen der weiche spongiöse Kern. Das ist die konzeptionelle Fracht dieses Kapitels. Eine vollständige Umbau-Rechnung, in der die Belastung das Material verändert, wäre eine zweite neue Idee — sie bekommt ein eigenes, noch folgendes Kapitel; hier erscheint sie nur als Ausblick.
16.2 Die 2D-Idealisierung, ehrlich
Ein echter Femur ist dreidimensional, gekrümmt und verdreht. Wir nehmen einen Längsschnitt durch die Frontalebene — Hüftkopf, Schenkelhals, Schaft — und vernetzen ihn mit allgemeinen P1-Dreiecken wie in Kapitel 10. Die runde Außenkontur nähern wir stückweise durch gerade Randkanten an; die Geometrie ist von Hand parametrisch gezeichnet, nicht aus einem CT-Bild importiert.
HinweisWas wir weglassen — und warum das Bild trotzdem stimmt
Wie beim „Rechteck statt I-Profil” aus Kapitel 1 nennen wir die Annahmen offen:
Ebener Spannungszustand. Wir rechnen mit derselben Materialtabelle wie in Kapitel 14 — dem ebenen Spannungszustand einer in die Tiefe freien Scheibe. Der Knochen ist das nicht ganz; aber das Buch lehrt bewusst nur diesen Zustand, und für ein qualitatives Lastbild ist die Vereinfachung brauchbar.
Ersatzscheibe fester Dicke. Der 2D-Schnitt steht für eine Scheibe der Dicke 0,02 m. Die volle Modellresultierende von 1,8 kN wird auf eine festgelegte Kopfkante und diese Ersatzdicke verteilt. Absolute Spannungen gelten deshalb nur für dieses Modell und sind keine klinische Aussage über einen echten Patienten.
Verrundeter Hals. Wo der Schenkelhals in den Schaft übergeht, ist ein echter Knochen verrundet — und wir runden dort ebenfalls. Eine scharfe einspringende Ecke erzeugte sonst eine Spannungssingularität: eine Spitze, die mit jedem feineren Netz weiterwächst und nie konvergiert. Die Verrundung macht die Halsspannung zu einer sinnvollen, konvergierenden Zahl. Die scharfe Ecke lassen wir später einmal bewusst auftreten — als Kerbwirkungs-Demo im Prüfteil.
16.3 Das ortsabhängige Material — die neue Idee
Wie kommt man vom Knochen zum E-Modul? Nicht in einem Schritt. Die ehrliche Messkette lautet:
Eine rohe Graustufe aus dem Röntgen- oder CT-Bild ist noch keine Dichte; erst eine kalibrierte Auswertung führt über die Aschedichte zur Scheindichte (Schileo u. a. 2008). Wir importieren kein CT, sondern verwenden ein synthetisches Scheindichtefeld, das nur die letzte Stufe nachbildet: dichter entlang des Hauptlastpfads (vom Kopf zum medialen Schaftfuß), dünner in den Randzonen.
Die letzte Stufe, von der Scheindichte zum E-Modul, ist das Herzstück. Für den spongiösen Knochen benutzen wir das Potenzgesetz von Morgan (Morgan u. a. 2003) (Dichte \(\rho\) in g/cm³, \(E\) in GPa):
\[
E \;=\; 6{,}95 \cdot \rho^{\,1{,}49}.
\]
Dieses Gesetz gilt nur im belegten spongiösen Bereich \(\rho \in [0{,}2;\,1{,}0]\) g/cm³. Wir extrapolieren es nicht bis zur kortikalen Dichte (rund 1,9 g/cm³) — die kortikale Schale bekommt separat ihren gemessenen Wert \(E = 17\) GPa. Morgan und Kollegen betonen außerdem, dass dieser Dichte-Steifigkeits-Zusammenhang vom anatomischen Ort abhängt: Die Formel ist ein didaktischer Arbeitskanon für den spongiösen Bereich, kein universelles Knochengesetz (Morgan u. a. 2003). Die folgende Zelle rechnet ein paar Belegwerte nach (Modul programme/gemeinsam/dichtefeld.py):
Man liest ab: Über den ganzen spongiösen Bereich reicht \(E\) von rund 0,6 bis knapp 7 GPa — und die kortikale Schale ist mit 17 GPa noch einmal gut doppelt so steif wie der dichteste Trabekelknochen. Der Knochen ist also nicht ein Material, sondern eine ganze Landkarte von Steifigkeiten. Dieselbe Kurve als Bild:
Abbildung 16.1: Das Morgan-Gesetz E = 6,95·ρ^1,49 im belegten spongiösen Bereich (0,2–1,0 g/cm³, durchgezogen). Die kortikale Schale (roter Punkt, fester Messwert 17 GPa bei ρ ≈ 1,9 g/cm³) wird bewusst NICHT über die verlängerte Kurve (gestrichelt) bestimmt: Morgans Daten decken den kortikalen Dichtebereich gar nicht ab, deshalb der separate Messwert.
HinweisEin historischer Vergleich: Carter & Hayes
Vor Morgan beschrieben Carter & Hayes (1977) die Steifigkeit aus Druckversuchen mit einer kubischen Dichteabhängigkeit, \(E \approx 3{,}79 \cdot
\rho^{3}\), zusätzlich schwach von der Dehnrate abhängig (Carter und Hayes 1977). Das Modul dichtefeld.py stellt beide Gesetze bereit — aber Carter & Hayes ist ein historischer Vergleich mit einer hier fest gewählten Dehnrate, kein universeller Materialschalter. Wir rechnen mit Morgan.
Die eigentliche Schlüsselgrafik dieses Kapitels ist diese Landkarte, über den vernetzten Femurschnitt gelegt:
Code
import os, syssys.path.insert(0, os.path.join("..", "programme", "gemeinsam"))sys.path.insert(0, os.path.join("..", "programme", "kap16"))import matplotlib.pyplot as pltfrom matplotlib.tri import Triangulationfrom dichtefeld import e_module_landkartefrom femur_netz import KNOTEN, DREIECKE, KORTIKAL, DICHTENe_module = e_module_landkarte(DICHTEN, KORTIKAL)xs = [k[0] for k in KNOTEN]ys = [k[1] for k in KNOTEN]tris = [list(d) for d in DREIECKE]e_gpa = [e /1e9for e in e_module]fig, ax = plt.subplots(figsize=(4.2, 6.2))tp = ax.tripcolor(Triangulation(xs, ys, tris), facecolors=e_gpa, cmap="viridis", edgecolors="0.6", linewidth=0.2)ax.set_aspect("equal")ax.set_xlabel("x (m)")ax.set_ylabel("y (m)")ax.set_title("E-Modul-Landkarte (GPa)")fig.colorbar(tp, ax=ax, shrink=0.7, label="E (GPa)")plt.tight_layout()plt.show()
Abbildung 16.2: Die Material-Landkarte des Femurschnitts: jedes Dreieck trägt seinen eigenen E-Modul. Die helle kortikale Randschale (17 GPa) umschließt den dunkleren spongiösen Kern, dessen Steifigkeit dem synthetischen Dichtefeld folgt — dichter (heller) entlang des Lastpfads vom Kopf zum medialen Schaftfuß.
16.4 Assemblieren mit variablem E
Was ändert sich am Programm? Erstaunlich wenig. Der Baustein aus Kapitel 14, baue_element_6x6, nahm den E-Modul schon immer pro Dreieck entgegen — wir haben ihn bisher nur für jedes Dreieck mit demselben Wert gefüttert. Jetzt geben wir jedem Dreieck seinen eigenen. Die einzige Änderung ist, dass die Materialtabelle je Dreieck aus dessen Dichte gebaut wird — eine Zeile.
Sehen wir es an einem Miniaturbeispiel: ein Streifen aus zwei Dreiecken, das linke steif, das rechte weich. Die Funktion assembliere_elastik_inhomogen nimmt eine Liste von E-Moduln, einen Wert je Dreieck:
Code
import os, syssys.path.insert(0, os.path.join("..", "programme", "gemeinsam"))from elastik2d import assembliere_elastik_inhomogen# Ein Rechteck-Streifen aus zwei Dreiecken.knoten = [(0.0, 0.0), (0.1, 0.0), (0.0, 0.05), (0.1, 0.05)]dreiecke = [(0, 1, 3), (0, 3, 2)]# Zwei verschiedene Materialien: links steif, rechts weich.e_module = [17.0e9, 2.0e9] # ein E-Modul JE Dreieckquerkontraktion =0.30tiefe =0.02matrix = assembliere_elastik_inhomogen(knoten, dreiecke, e_module, querkontraktion, tiefe)# Das steifste Diagonalelement heraussuchen (reines Python, keine Kurzform).groesstes =0.0for i inrange(len(matrix)):if matrix[i][i] > groesstes: groesstes = matrix[i][i]print(f"Systemmatrix: {len(matrix)} x {len(matrix)} (4 Knoten x 2 Richtungen)")print(f"steifstes Diagonalelement: {groesstes /1e6:8.1f} MN/m")
Systemmatrix: 8 x 8 (4 Knoten x 2 Richtungen)
steifstes Diagonalelement: 406.3 MN/m
Der Assemblierungs-Loop ist Wort für Wort der aus Kapitel 10 und Kapitel 14 — für jedes Dreieck den 6×6-Steckbrief bauen und an die Freiheitsgrade seiner drei Knoten addieren. Neu ist allein, dass der Steckbrief mit e_module[nummer] statt mit einem globalen e_modul gebaut wird. Das Kapitelprogramm ruft danach genau diese geprüfte Funktion über alle 427 Dreiecke des Femurs auf, statt den Loop zu wiederholen.
WichtigAuch das Spannungslesen braucht das lokale E
Es gibt eine zweite Stelle, an der das ortsabhängige Material zählt: das Zurückrechnen der Spannungen. Eine steife kortikale Schale trägt bei gleicher Dehnung mehr Spannung als der weiche Kern. Würde man die Spannungen mit einem globalen \(E\) auslesen, stimmte das Verhältnis der Felder nicht. Deshalb gibt es neben assembliere_elastik_inhomogen auch elementspannungen_inhomogen, die jedem Dreieck seine eigene Materialtabelle gibt. Beide Funktionen sind so gebaut, dass sie bei überall gleichem \(E\) Zeile für Zeile dasselbe liefern wie die vertrauten Funktionen aus Kapitel 14.
16.5 Lasten, Lager — und die Rechnung
Die Hüftkontaktkraft von 1,8 kN wirkt nach unten auf eine kurze Kante der Kopfoberseite. Wie in Kapitel 13 verteilen wir eine Kantenlast in äquivalente Knotenkräfte: Jede Randkante gibt ihren Anteil hälftig an ihre beiden Endknoten. Der Schaftfuß bei \(y = 0\) wird in beide Richtungen festgehalten.
HinweisDer schwimmende Knochen
Wie in Kapitel 8 und Kapitel 14 gilt: Ohne Auflager schwimmt auch der Knochen davon — die Systemmatrix hätte Starrkörpermoden und wäre nicht lösbar. Das künstliche Auflager am Schaftfuß sperrt diese Moden. Es steht für „hier geht der Knochen weiter, den wir nicht mitmodellieren”, nicht für ein echtes Gelenk.
Jetzt läuft die vertraute Kette. Die folgende Zelle assembliert das ganze Femurnetz mit seiner Material-Landkarte, bringt die Hüftlast auf, hält den Schaftfuß fest, löst — und liest die Kennzahlen aus. Sie benutzt nur geprüfte Funktionen aus elastik2d.py:
Code
import os, syssys.path.insert(0, os.path.join("..", "programme", "gemeinsam"))sys.path.insert(0, os.path.join("..", "programme", "kap16"))from elastik2d import (assembliere_elastik_inhomogen, elementspannungen_inhomogen, loese_elastik)from dichtefeld import e_module_landkartefrom femur_netz import KNOTEN, DREIECKE, KORTIKAL, DICHTENQUERKONTRAKTION =0.30TIEFE =0.02HUEFTLAST =1800.0LAST_Y =0.155# 1) Material-Landkarte: ein E-Modul je Dreieck.e_module = e_module_landkarte(DICHTEN, KORTIKAL)# 2) Systemmatrix mit variablem E.matrix = assembliere_elastik_inhomogen(KNOTEN, DREIECKE, e_module, QUERKONTRAKTION, TIEFE)# 3) Hueftlast als aequivalente Knotenkraefte auf die Kopfkante (nach unten).# Erst die Knoten auf der Kopfoberkante sammeln, dann nach x sortieren# (ueber Paare (x, Knotennummer) -- reines Python, ohne Lambda).kanten_paare = []for i inrange(len(KNOTEN)): x = KNOTEN[i][0] y = KNOTEN[i][1]ifabs(y - LAST_Y) <1e-6and-0.010-1e-9<= x <=0.030+1e-9: kanten_paare.append((x, i))kanten_paare.sort()auf_kante = []for paar in kanten_paare: auf_kante.append(paar[1])traktion = HUEFTLAST / ((0.030- (-0.010)) * TIEFE)kraefte = [0.0] * (2*len(KNOTEN))for stelle inrange(len(auf_kante) -1): links = auf_kante[stelle] rechts = auf_kante[stelle +1] halbe = traktion *abs(KNOTEN[rechts][0] - KNOTEN[links][0]) * TIEFE /2.0 kraefte[2* links +1] = kraefte[2* links +1] - halbe kraefte[2* rechts +1] = kraefte[2* rechts +1] - halbe# 4) Schaftfuss festhalten und loesen.lager = {}for i inrange(len(KNOTEN)):if KNOTEN[i][1] <1e-9: lager[2* i] =0.0 lager[2* i +1] =0.0verschiebungen = loese_elastik(matrix, kraefte, lager)# 5) Spannungen mit dem jeweiligen Element-E (rein mechanisch: keine Waerme).null_dt = [0.0] *len(DREIECKE)spannungen = elementspannungen_inhomogen(KNOTEN, DREIECKE, verschiebungen, e_module, QUERKONTRAKTION, 0.0, null_dt)# Kennzahlen: Kopfsenkung und flaechengewichtete Halsspannung.tiefste =0.0for i inrange(len(KNOTEN)):if verschiebungen[2* i +1] < tiefste: tiefste = verschiebungen[2* i +1]kopfsenkung_mm =-tiefste *1000.0def von_mises(spannung): sxx, syy, txy = spannungreturn (sxx * sxx - sxx * syy + syy * syy +3.0* txy * txy) **0.5summe =0.0flaeche_hals =0.0for nummer inrange(len(DREIECKE)): d = DREIECKE[nummer] sy = (KNOTEN[d[0]][1] + KNOTEN[d[1]][1] + KNOTEN[d[2]][1]) /3.0if0.082< sy <0.106: (x0, y0), (x1, y1), (x2, y2) = KNOTEN[d[0]], KNOTEN[d[1]], KNOTEN[d[2]] flaeche =0.5*abs((x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0)) summe = summe + von_mises(spannungen[nummer]) * flaeche flaeche_hals = flaeche_hals + flaechehalsspannung_mpa = summe / flaeche_hals /1e6print(f"Kopfsenkung (grobes Lehrnetz) : {kopfsenkung_mm:.3f} mm")print(f"Halsspannung, flaechengewichtet : {halsspannung_mpa:.2f} MPa")
Zwei Zahlen zum Merken: Der Hüftkopf senkt sich um rund ein Zehntel Millimeter, und im Schenkelhals herrscht eine gemittelte Von-Mises-Vergleichsspannung von gut 5 MPa. Das grobe Lehrnetz rechnet der reine-Python-Löser in Sekunden.
WichtigGrobe Netze sind zu steif — der Fall für Kapitel 12
Die hier gezeigte Kopfsenkung (rund 0,12 mm) ist noch nicht die ganze Wahrheit. Wie in Kapitel 12 und Kapitel 14 sind grobe P1-Netze zu steif: Sie unterschätzen die Verformung und nähern sich der richtigen Lösung von unten. Verfeinert man das Netz systematisch, konvergieren Kopfsenkung und gemittelte Halsspannung gegen die validierten Referenzwerte 0,15 mm und 5,0 MPa (nachgerechnet mit dem professionellen Löser scikit-fem auf immer feineren Netzen — Protokoll in validierung/kap16_knochen/). Die grobe Zahl ist nicht falsch, nur noch nicht konvergiert.
Das Ergebnisbild zeigt links die Material-Landkarte noch einmal, in der Mitte das überhöht verformte Netz und rechts das Von-Mises-Feld:
Code
import matplotlib.pyplot as pltfrom matplotlib.tri import Triangulatione_gpa = [e /1e9for e in e_module]vm_mpa = [von_mises(s) /1e6for s in spannungen]ueberhoehung =20.0vx = [KNOTEN[i][0] + ueberhoehung * verschiebungen[2* i] for i inrange(len(KNOTEN))]vy = [KNOTEN[i][1] + ueberhoehung * verschiebungen[2* i +1] for i inrange(len(KNOTEN))]tris = [list(d) for d in DREIECKE]xs = [k[0] for k in KNOTEN]ys = [k[1] for k in KNOTEN]fig, achsen = plt.subplots(1, 3, figsize=(9.6, 6.0))tpe = achsen[0].tripcolor(Triangulation(xs, ys, tris), facecolors=e_gpa, cmap="viridis", edgecolors="face")achsen[0].set_title("E-Landkarte (GPa)")fig.colorbar(tpe, ax=achsen[0], shrink=0.7)achsen[1].triplot(xs, ys, tris, color="0.8", linewidth=0.4)achsen[1].triplot(vx, vy, tris, color="C3", linewidth=0.4)achsen[1].set_title("verformt (×20 überhöht)")tp = achsen[2].tripcolor(Triangulation(vx, vy, tris), facecolors=vm_mpa, cmap="inferno", edgecolors="face")achsen[2].set_title("Von-Mises-Spannung (MPa)")fig.colorbar(tp, ax=achsen[2], shrink=0.7)for ax in achsen: ax.set_aspect("equal") ax.set_xlabel("x (m)")achsen[0].set_ylabel("y (m)")plt.tight_layout()plt.show()
Abbildung 16.3: Ergebnis des Femurmodells, drei Bilder. Links die E-Modul-Landkarte (kortikale Schale 17 GPa, spongiöser Kern nach Dichtefeld). Mitte: das überhöht (×20) verformte Netz über der Ausgangslage — der Kopf senkt sich, der Hals biegt sich leicht. Rechts die Von-Mises-Vergleichsspannung — hoch am Schenkelhals, wo der Querschnitt am schmalsten ist, und in der lasttragenden kortikalen Schale. Von Mises ist hier nur ein Farbwert, kein Versagenskriterium.
Lesen ohne Übertreiben. Das Bild zeigt hohe Spannung am Schenkelhals und in der kortikalen Schale. Verlockend wäre der Satz „die Schale trägt alles” — aber hohe Spannung allein beweist noch keinen Lastanteil. Von Mises fasst nur die drei Spannungskomponenten zu einem skalaren Farbwert zusammen; es ist hier ein Darstellungswert, kein Versagenskriterium und keine Aussage darüber, ab wann ein Knochen bricht.
WichtigEchte Spitze oder Kerbwirkung? — die Kap.-12-Frage
Am Übergang vom Hals zum Kopf leuchten im Von-Mises-Bild kleine Spitzen auf. Sind die echt? Die Prüfung ist die aus Kapitel 12: Netz verfeinern und nachschauen. Für die verrundete Kontur konvergieren Kopfsenkung und gemittelte Halsspannung — sie sind belastbar. Für eine scharfe, unverrundete Ecke dagegen wächst die Spitzenspannung mit jedem feineren Netz weiter (in der Nachrechnung von rund 9 über 12 auf über 13 MPa und weiter): Das ist eine mathematische Singularität, kein physikalischer Messwert. Genau darum haben wir den Hals verrundet — damit die Halsspannung eine Zahl ist, der man trauen kann, und nicht ein Artefakt der Netzweite.
Der Lastaufbau, Bild für Bild. Weil das Modell linear ist, müssen wir nicht für jede Laststufe neu rechnen: Wächst die Hüftkraft von null auf 1,8 kN, so wachsen Verformung und Von-Mises-Feld proportional mit. Die folgende Darstellung skaliert die eine berechnete Lösung — der Kopf sinkt tiefer, das Spannungsfeld leuchtet stärker auf, aber das Muster bleibt sich gleich. Genau das bedeutet Linearität.
Code
import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.tri import Triangulationfrom matplotlib.animation import FuncAnimationfrom matplotlib.cm import ScalarMappablefrom matplotlib.colors import Normalizefrom IPython.display import HTMLbasis_vm = [von_mises(s) /1e6for s in spannungen]vmax =max(basis_vm)xs = [k[0] for k in KNOTEN]ys = [k[1] for k in KNOTEN]tris = [list(d) for d in DREIECKE]ueberhoehung =20.0faktoren = [i /23.0for i inrange(24)]norm = Normalize(vmin=0.0, vmax=vmax)fig, ax = plt.subplots(figsize=(4.2, 6.2))sm = ScalarMappable(norm=norm, cmap="inferno")sm.set_array([])fig.colorbar(sm, ax=ax, shrink=0.7, label="Von Mises (MPa)")def zeichne(k): ax.clear() f = faktoren[k] vx = [xs[i] + ueberhoehung * f * verschiebungen[2* i] for i inrange(len(KNOTEN))] vy = [ys[i] + ueberhoehung * f * verschiebungen[2* i +1] for i inrange(len(KNOTEN))] vm = [f * wert for wert in basis_vm] ax.tripcolor(Triangulation(vx, vy, tris), facecolors=vm, cmap="inferno", norm=norm, edgecolors="face") ax.set_xlim(-0.02, 0.045) ax.set_ylim(-0.01, 0.165) ax.set_aspect("equal") ax.set_title("Hüftkraft %5.0f N"% (f *1800.0)) ax.axis("off")return []ani = FuncAnimation(fig, zeichne, frames=len(faktoren), interval=90, blit=False)plt.close(fig)HTML(ani.to_jshtml())
Abbildung 16.4: Der Lastaufbau am Femur: Die Hüftkontaktkraft wächst von 0 auf 1,8 kN, das überhöht gezeichnete Netz senkt sich und das Von-Mises-Feld baut sich proportional auf. Weil das Modell linear ist, ist jedes Standbild die auf die jeweilige Kraft skalierte Endlösung — kein neuer Rechenlauf, nur ein Faktor.
Was man sieht: Bei kleiner Kraft ist alles dunkel und fast unverformt; mit wachsender Last senkt sich der Kopf und der Schenkelhals leuchtet zuerst auf. Nichts an der Verteilung ändert sich — sie ist von Anfang an da, nur schwächer. Ein nichtlineares Material (oder große Verformungen) würde hier ein anderes Muster bei starker Last zeigen; unser linear-elastischer Knochen tut das bewusst nicht.
16.6 Probier es aus: die Knochenwerkstatt
Zeit für dein eigenes Experiment. Die folgende interaktive Einheit rechnet ein grobes Femurnetz im Browser. Du kannst zwei Dinge einstellen: den Faktor der Hüftkraft (in Vielfachen des Körpergewichts) und die Steifigkeit des spongiösen Kerns über eine repräsentative Dichte.
WichtigVorhersage-Punkt
Bevor du rechnest: Wenn die Spongiosa dichter und damit steifer wird — wird die Kopfsenkung dann größer oder kleiner? Leg dich fest, dann probier es aus.
Abbildung 16.5: Vorgerenderte Fassung der Knochenwerkstatt in den Ausgangseinstellungen (Hüftkraft 2,5·BW, Spongiosa-Dichte 0,5 g/cm³): E-Landkarte, überhöht verformtes Netz und Von-Mises-Feld auf demselben groben 95-Knoten-Netz wie das interaktive Stück. Im Browser wird dieses Bild durch dein eigenes Ergebnis ersetzt.
Die Antwort auf den Vorhersage-Punkt: Dichtere Spongiosa ist steifer, also senkt sich der Kopf weniger. Weil das Modell linear ist, verdoppelt eine doppelte Hüftkraft dagegen jede Verschiebung und jede Spannung exakt — das kannst du direkt am Kraftfaktor ablesen.
16.7 Ausblick: Warum Knochen sich umbauen
Wir haben aus einer festen Dichtekarte eine E-Modul-Karte gemacht und den Knochen einmal durchgerechnet. Aber Knochen ist lebendig: Über Wochen und Monate baut er dort Material auf, wo er stark belastet wird, und dort ab, wo er es kaum ist. Diese Beobachtung trägt einen Namen — das Wolffsche Gesetz — und einen modernen biologischen Rahmen, Frosts Mechanostat(Frost 2003).
Moderne FE-Remodelingmodelle gießen das in eine Rückkopplung: Sie vergleichen einen mechanischen Reiz (etwa die Formänderungsenergiedichte) mit einem Zielwert — oberhalb wird Material aufgebaut, unterhalb abgebaut, in einer Totzone dazwischen bleibt alles, wie es ist (Huiskes u. a. 1987). Das ist genau die Form einer Rückkopplung, deren Tür 1 aufgestoßen hat: Diesmal reagiert nicht die Temperatur auf die Wärme, sondern das Material auf die Rechnung.
Anders als der Einweg-Pfeil des Trägers schließt sich dieser Kreis:
Code
import numpy as npimport matplotlib.pyplot as pltnamen = ["Dichte ρ", "E-Modul\n(Morgan)", "FEM-\nRechnung","mechan.\nReiz", "neue\nDichte"]winkel = np.linspace(90, 90-360, len(namen), endpoint=False) * np.pi /180.0rx, ry =1.0, 1.0fig, ax = plt.subplots(figsize=(6.2, 5.2))px = rx * np.cos(winkel)py = ry * np.sin(winkel)for i inrange(len(namen)): ax.text(px[i], py[i], namen[i], ha="center", va="center", fontsize=10, bbox=dict(boxstyle="round,pad=0.4", fc="#dfe8f0", ec="black"))for i inrange(len(namen)): j = (i +1) %len(namen) ax.annotate("", xy=(0.72* px[j], 0.72* py[j]), xytext=(0.72* px[i], 0.72* py[i]), arrowprops=dict(arrowstyle="-|>", color="tab:red", lw=1.8, connectionstyle="arc3,rad=0.2"))ax.set_xlim(-1.7, 1.7)ax.set_ylim(-1.6, 1.6)ax.set_aspect("equal")ax.axis("off")plt.tight_layout()plt.show()
Abbildung 16.6: Die Remodeling-Rückkopplung als geschlossener Kreis: Aus der Dichte wird über Morgan ein E-Modul, damit rechnet die FEM Verschiebungen und Spannungen, daraus ein mechanischer Reiz, und der verändert die Dichte für die nächste Runde. Genau dieser geschlossene Kreis unterscheidet das Remodeling von der Einweg-Kopplung des Trägers — gerechnet wird er erst im folgenden Kapitel.
In diesem Kapitel rechnen wir das nicht. Eine Rückkopplung will sorgfältig gezähmt werden: Ohne Begrenzung, Totzone und räumliche Mittelung erzeugt eine naive Update-Regel bunte Schachbrettmuster, die wie Trabekel aussehen, aber nur numerische Artefakte sind. Wie man die Schleife baut und ihr misstraut, ist die Aufgabe des folgenden Kapitels über das Remodeling.
TippMerkkasten
Ein tragender Knochen ist für kleine, kurze Lasten dasselbe FEM-Problem wie ein Träger — dieselbe Maschine (elastik2d.py), rein mechanisch.
Die Hüftkontaktkraft (~2,5·BW ≈ 1,8 kN) ist eine Gelenkresultierende und enthält die Muskelwirkung bereits — kein zweiter Muskel-Lastfall.
Das Material ist eine Landkarte: kortikale Schale 17 GPa, spongiöser Kern aus dem Dichtefeld über Morgan (\(E = 6{,}95\,\rho^{1{,}49}\), nur \(0{,}2\)–\(1{,}0\) g/cm³).
Zwei Funktionen brauchen das lokale E: Assemblieren und Spannungslesen.
Von Mises ist hier ein Farbwert, kein Versagenskriterium; scharfe Ecken erzeugen Kerbsingularitäten — verrunden, dann konvergiert die Halsspannung.
Roter Faden
Zurück: die 2D-Elastik-Kette \(u \to \varepsilon \to \sigma\) (Kapitel 14), Auflager und Starrkörpermoden (Kapitel 8, Kapitel 14), die Knotenkraft (Kapitel 13), die Fachwerk-Analogie (Kapitel 3), die Prüfrituale (Kapitel 12) und die Rückkopplungs-Tür (1). — Vor: das folgende Kapitel schließt die Remodeling-Schleife; die Übertragen-Übung unten weist über den Knochen hinaus in Zahnmedizin und Botanik.
Übungen
Ü 16.1 (Verstehen). Im E-Landkarten-Bild (Abbildung 16.2) sind zwei Zonen zu sehen. Welche ist die kortikale Schale, welche der spongiöse Kern? Und warum müssen E-Feld, Dehnungsfeld und Spannungsfeld nicht dasselbe Bild zeigen?
HinweisMusterlösung zu Ü 16.1
Die helle, gleichmäßig steife Randzone (17 GPa) ist die Kompakta; der dunklere, fleckige Innenbereich ist die Spongiosa (0,6–7 GPa je nach Dichte). Die drei Felder unterscheiden sich, weil sie verschiedene Dinge zeigen: \(E\) ist eine reine Materialeigenschaft (wo ist der Knochen steif?), die Dehnung misst, wie stark sich das Material verformt, und die Spannung ist deren Produkt mit \(E\) — eine steife Zone kann bei kleiner Dehnung hohe Spannung tragen. Deshalb sitzt die höchste Spannung nicht dort, wo \(E\) am größten ist, sondern wo Last und Steifigkeit zusammentreffen: am schmalen Schenkelhals und in der Schale.
Ü 16.2 (Verändern). Verdopple in der Knochenwerkstatt den Kraftfaktor von 2 auf 4. Verdoppelt sich die Kopfsenkung? Stelle danach die Spongiosa-Dichte von 0,3 auf 0,9 g/cm³ und vergleiche.
HinweisMusterlösung zu Ü 16.2
Ja — die Kopfsenkung verdoppelt sich exakt, denn das Modell ist linear: Alle Verschiebungen und Spannungen wachsen proportional zur Last. Eine dichtere Spongiosa (0,9 statt 0,3 g/cm³) hebt deren E-Modul von rund 1,2 auf gut 5,9 GPa; der Kern wird steifer, die Kopfsenkung kleiner. Nichtlinear wäre die Dichte nur, wenn sie über Morgans Potenz \(\rho^{1{,}49}\) ins E einginge — was sie tut, aber die Rechnung selbst bleibt bei fester Landkarte linear.
Ü 16.3 (Übertragen). Suche dir ein zweites tragendes Bauteil aus der Biologie — ein Zahnimplantat im Kieferknochen oder einen pflanzlichen Stängel als Rohr-Träger. Benenne (ohne zu rechnen): Was ist das Netz, was die Last, was das Auflager, und wo säße das ortsabhängige Material?
HinweisMusterlösung zu Ü 16.3
Zahnimplantat: Das Netz ist ein Schnitt durch Implantat und umgebenden Kieferknochen; die Last ist die Kaukraft auf die Krone; das Auflager ist der Rand des mitmodellierten Knochenblocks. Ortsabhängiges Material: hartes Titan im Implantat, kortikale Schicht und weiche Spongiosa im Kiefer — dieselbe Landkarten-Idee. Pflanzenstängel: Netz ist ein Längs- oder Querschnitt; Last ist Wind oder Eigengewicht der Blüte; Auflager ist der Wurzelansatz. Das Material ist ortsabhängig zwischen fester Randfaser (Sklerenchym) und weichem Mark — ein Rohr-Träger wie das I-Profil aus Kapitel 1.
Das Kleingedruckte
Die Grenzen dieses Modells, offen benannt: Es ist 2D, nicht 3D — ein echter Femur trägt räumlich. Es ist isotrop, obwohl Knochen faserig und richtungsabhängig ist. Es ist eine statische Ersatzrechnung eines einzigen Lastmoments, kein Wachstum, keine Heilung, kein Umbau. Von Mises ist ein Darstellungswert, kein Bruchkriterium; wir treffen keine klinische Aussage über einen realen Knochen. Und die Geometrie ist von Hand gezeichnet, nicht aus einem CT importiert. Jede dieser Grenzen ist eine bewusste Vereinfachung — genau die Ehrlichkeit, die Kapitel 12 verlangt. Was das Modell dennoch zeigt, ist echt: dass dieselbe Methode, die den Träger im Feuer durchrechnet, ohne eine neue Idee auch einen Knochen trägt.
Bergmann, Georg, G. Deuretzbacher, M. Heller, u. a. 2001. „Hip contact forces and gait patterns from routine activities“. Journal of Biomechanics 34 (7): 859–71. https://doi.org/10.1016/S0021-9290(01)00040-9.
Carter, Dennis R., und Wilson C. Hayes. 1977. „The compressive behavior of bone as a two-phase porous structure“. The Journal of Bone and Joint Surgery 59 (7): 954–62. https://doi.org/10.2106/00004623-197759070-00021.
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.
Morgan, Elise F., Harun H. Bayraktar, und Tony M. Keaveny. 2003. „Trabecular bone modulus-density relationships depend on anatomic site“. Journal of Biomechanics 36 (7): 897–904. https://doi.org/10.1016/S0021-9290(03)00071-X.
Schileo, Enrico, Enrico Dall’Ara, Fulvia Taddei, u. a. 2008. „An accurate estimation of bone density improves the accuracy of subject-specific finite element models“. Journal of Biomechanics 41 (11): 2483–91. https://doi.org/10.1016/j.jbiomech.2008.05.017.
