11  Von Hand bis zum Bild: das Handnetz vollständig gerechnet

Am Ende von Kapitel 10 stand ein zusammengebautes Netz: dreizehn Dreiecke über der Ferse, Element für Element zu einer großen Steifigkeitsmatrix verbunden. Aber ein Netz, das weder belastet noch gelagert ist, gibt noch kein Bild; es schwebt. Dieses Kapitel holt nach, was der Name Handnetz verspricht: Es rechnet das Netz von Hand durch, vom Zusammenbau der Matrix über die Last und das Lager bis zum ersten Feldbild. Aus der Assemblierung, dem Lastvektor eines Schrittmoments und dem Talus-Lager entsteht das erste Spannungsbild, dessen jede Zahl du selbst nachrechnen kannst.

Lernziele

Nach diesem Kapitel kannst du

  1. eine Systemmatrix als Zusammenbau aller Elementmatrizen lesen und ihrem Besetzungsmuster ansehen, wie das Netz gebaut ist,
  2. eine über den Rand verteilte Kraft in Knotenkräfte übersetzen und einem fertigen Lastvektor ansehen, welche Kraft wo angreift,
  3. ein festes Lager in Matrix und Lastvektor einbauen — einmal in der vollen Matrix, einmal durch Streichen — und begründen, warum beides dasselbe tut,
  4. aus einer gelösten Verschiebung ein Spannungsfeld gewinnen und sagen, was die Von-Mises-Farbe zeigt und was sie verschweigt,
  5. die Kraft ausrechnen, mit der ein Lager hält, und mit den aufgebrachten Kräften ins Gleichgewicht setzen,
  6. beurteilen, welche Aussagen ein grobes Netz trägt und welche nicht.

11.1 Der Fahrplan: von der Matrix zum Bild

Vier Schritte trennen das schwebende Netz vom fertigen Spannungsbild, und es sind genau die vier, die jede finite Elemente-Rechnung geht — nur hier so klein, dass jeder Schritt einzeln sichtbar wird. Zuerst die Systemmatrix: die dreizehn Elementmatrizen an ihre Plätze einsortiert. Dann der Lastvektor: die drei Kräfte des Schrittmoments auf die Randknoten verteilt. Dann das Lager: die Talus-Knoten festgehalten, damit die Ferse nicht davonschwebt. Und schließlich die Lösung: das lineare Gleichungssystem gelöst, die Verschiebung in ein Spannungsbild übersetzt.

Damit das Kapitel für sich läuft, stehen die Bausteine aus Kapitel 10 — die Elementsteifigkeitsmatrix eines Dreiecks und die Assemblierung — hier noch einmal gefaltet. Neu ist alles, was danach kommt. Wer den folgenden Codeblock überspringt, verliert keinen Gedanken: er wiederholt nur die Rechenbausteine des vorigen Kapitels und legt die festen Zahlen dieses Kapitels fest. Material, Scheibendicke und Körpergewicht holt er dabei aus dem Zahlenanhang (Anhang A) statt sie hier noch einmal hinzuschreiben; Handnetz, Randgruppen und die fünf Momente stehen darunter mit ihrer Herkunft im Kommentar.

Das Rechenwerk dieses Kapitels — die Bausteine aus Kapitel 10 zum Nachlesen aufklappen 
import math
import sys

sys.path.insert(0, "programme/gemeinsam")
import kanon

# --- Material aus dem Kanon -------------------------------------------------
E_KORTIKAL = kanon.wert_von("E_KORTIKAL")   # Pa, kortikaler Knochen, homogen
NU_KNOCHEN = kanon.wert_von("NU")           # Querdehnzahl
TIEFE = kanon.wert_von("DICKE")             # m, Ersatzdicke der 2D-Scheibe
KG = kanon.wert_von("KOERPERGEWICHT")       # N, Körpergewicht der Referenzperson

# --- Die fünf Schrittmomente (identisch mit denen der Referenzrechnung) -----
# (Name, Boden, Achilles, Faszie) in Vielfachen von KG.
MOMENTE = [
    ("Fersenauftritt",      1.45, 0.30, 0.10),
    ("Belastungsaufbau",    2.00, 1.50, 0.30),
    ("mittlere Standphase", 2.50, 2.00, 0.60),
    ("Fersenabloesung",     1.20, 3.50, 1.20),
    ("Abstoss",             0.30, 2.80, 2.00),
]
RICHTUNG_BODEN = (0.0, 1.0)          # Druck nach oben auf die plantare Flaeche
RICHTUNG_ACHILLES = (-0.30, 0.954)   # Zug nach hinten-oben (Achillessehne)
RICHTUNG_FASZIE = (0.985, -0.174)    # Zug nach vorn-unten (Plantarfaszie)

# --- Das Handnetz (identisch mit dem der Referenzrechnung) ------------------
# Doppelfaecher, 13 Knoten, 13 Dreiecke; innere Knoten 8 (links), 11 (rechts).
HANDNETZ_KNOTEN = [
    (0.006, 0.0000), (0.033, 0.0134), (0.052, 0.006), (0.066, 0.018),
    (0.050, 0.033), (0.030, 0.032), (0.012, 0.032), (0.001, 0.012),
    (0.016, 0.015), (0.0218, 0.0000), (0.0259, 0.0071), (0.036, 0.024),
    (0.004, 0.022)]
HANDNETZ_DREIECKE = [
    (7, 0, 8), (0, 9, 8), (9, 10, 8), (10, 1, 8),
    (5, 6, 8), (6, 12, 8), (12, 7, 8),
    (1, 2, 11), (2, 3, 11), (3, 4, 11), (4, 5, 11),
    (1, 11, 8), (5, 8, 11)]

# --- Randgruppen aus dem Kanon (identisch mit der Referenzrechnung) ---------
BODEN_KANTEN = [(0, 7), (0, 9), (9, 10)]     # plantare Bodenkontaktfacetten
ENTHESE_KANTEN = [(9, 10)]                   # Enthese am Sporn-Ort (Faszie)
ACHILLES_KANTEN = [(12, 7)]                  # Achillessehnen-Ansatz
TALUS_LAGERKNOTEN = [4, 5, 6]                # Talus-Gelenkflaeche (festes Lager)

def zeltsteigungen(ecken):
    """Die konstanten Steigungen b, c der drei Zeltdach-Ansaetze und die Flaeche."""
    (xa, ya), (xb, yb), (xc, yc) = ecken[0], ecken[1], ecken[2]
    zwei = xa * (yb - yc) + xb * (yc - ya) + xc * (ya - yb)
    if zwei == 0.0:
        raise ValueError("Dreieck entartet (Flaeche null)")
    b = [(yb - yc) / zwei, (yc - ya) / zwei, (ya - yb) / zwei]
    c = [(xc - xb) / zwei, (xa - xc) / zwei, (xb - xa) / zwei]
    return b, c, abs(zwei) / 2.0

def materialtabelle(e_modul, nu):
    """Die 3x3-Materialtabelle des ebenen Spannungszustands (Hookesches Gesetz)."""
    faktor = e_modul / (1.0 - nu * nu)
    schub = (1.0 - nu) / 2.0
    return [[faktor, faktor * nu, 0.0],
            [faktor * nu, faktor, 0.0],
            [0.0, 0.0, faktor * schub]]

def b_tabelle(ecken):
    b, c, _ = zeltsteigungen(ecken)
    tab = [[0.0] * 6, [0.0] * 6, [0.0] * 6]
    for e in range(3):
        tab[0][2 * e] = b[e]
        tab[1][2 * e + 1] = c[e]
        tab[2][2 * e] = c[e]
        tab[2][2 * e + 1] = b[e]
    return tab

def elementmatrix(ecken, e_modul, nu, tiefe):
    """Die 6x6-Elementsteifigkeitsmatrix eines Dreiecks (Kapitel 10)."""
    B = b_tabelle(ecken)
    C = materialtabelle(e_modul, nu)
    _, _, flaeche = zeltsteigungen(ecken)
    cb = [[0.0] * 6 for _ in range(3)]
    for z in range(3):
        for s in range(6):
            summe = 0.0
            for k in range(3):
                summe = summe + C[z][k] * B[k][s]
            cb[z][s] = summe
    Ke = [[0.0] * 6 for _ in range(6)]
    for i in range(6):
        for j in range(6):
            summe = 0.0
            for k in range(3):
                summe = summe + B[k][i] * cb[k][j]
            Ke[i][j] = summe * tiefe * flaeche
    return Ke

def leeres_system(dof):
    return [[0.0] * dof for _ in range(dof)]

def assembliere(knoten, dreiecke, e_modul, nu, tiefe):
    """Baut aus den Elementmatrizen die grosse Systemmatrix (Kapitel 10)."""
    K = leeres_system(2 * len(knoten))
    for dr in dreiecke:
        ecken = [knoten[dr[0]], knoten[dr[1]], knoten[dr[2]]]
        Ke = elementmatrix(ecken, e_modul, nu, tiefe)
        dofs = []
        for e in range(3):
            dofs.append(2 * dr[e])
            dofs.append(2 * dr[e] + 1)
        for a in range(6):
            for bb in range(6):
                K[dofs[a]][dofs[bb]] = K[dofs[a]][dofs[bb]] + Ke[a][bb]
    return K

def gauss_loese(A, rhs):
    """Loest ein lineares Gleichungssystem mit Gauss-Elimination (Kapitel 5)."""
    n = len(rhs)
    a = [zeile[:] for zeile in A]
    b = rhs[:]
    for s in range(n):
        p = s
        for r in range(s + 1, n):
            if abs(a[r][s]) > abs(a[p][s]):
                p = r
        if a[p][s] == 0.0:
            raise ValueError("System nicht loesbar (Pivot null)")
        if p != s:
            a[s], a[p] = a[p], a[s]
            b[s], b[p] = b[p], b[s]
        for r in range(s + 1, n):
            faktor = a[r][s] / a[s][s]
            if faktor != 0.0:
                for c in range(s, n):
                    a[r][c] = a[r][c] - faktor * a[s][c]
                b[r] = b[r] - faktor * b[s]
    x = [0.0] * n
    for r in range(n - 1, -1, -1):
        summe = b[r]
        for c in range(r + 1, n):
            summe = summe - a[r][c] * x[c]
        x[r] = summe / a[r][r]
    return x

def feste_dofs(knoten_liste):
    """Die Freiheitsgrade (x und y) einer Liste festgehaltener Knoten."""
    feste = []
    for k in knoten_liste:
        feste.append(2 * k)
        feste.append(2 * k + 1)
    return feste

def komma(zahl, stellen=1):
    """Eine Zahl mit deutschem Dezimalkomma setzen (fuer den Fliesstext)."""
    return ("%.*f" % (stellen, zahl)).replace(".", ",")

11.2 Das Netz, auf dem gerechnet wird

Bevor die erste Matrix entsteht, lohnt ein Blick auf die Ausgangslage. Das Handnetz ist dasselbe wie in Kapitel 10 — dreizehn nummerierte Knoten, dreizehn Dreiecke, mit denselben Nummern wie dort. Neu ist allein, dass jeder Knoten jetzt eine von drei Rollen bekommt: Er trägt eine Kraft, er ist festgehalten, oder er ist frei und bewegt sich nur, weil seine Nachbarn ihn ziehen. Welche Rolle ein Knoten hat, entscheidet nicht das Kapitel, sondern die anatomische Einteilung des Randes: Sie legt fest, welche Randkanten Bodenkontakt sind, welche die Faszien-Enthese tragen, an welcher die Achillessehne zieht und welche auf der Talus-Gelenkfläche liegen. Abbildung 11.1 zeichnet das Netz mit genau dieser Einteilung.

Code 
import numpy as np
import matplotlib.pyplot as plt

kn = np.array(HANDNETZ_KNOTEN) * 1000.0        # m -> mm
RANDWEG = [0, 9, 10, 1, 2, 3, 4, 5, 6, 12, 7]  # geschlossener Randpolygonzug
TALUS_KANTEN = [(4, 5), (5, 6)]

fig, ax = plt.subplots(figsize=(6.8, 4.6))
for dr in HANDNETZ_DREIECKE:
    ax.fill(kn[list(dr), 0], kn[list(dr), 1], color="#f6efe4",
            ec="#bfa887", lw=1.0, zorder=1)

def kante(i, j, farbe, breite, stil="-", beschriftung=None):
    ax.plot([kn[i, 0], kn[j, 0]], [kn[i, 1], kn[j, 1]], stil, color=farbe,
            lw=breite, solid_capstyle="round", zorder=3, label=beschriftung)

# unklassifizierte Randkanten zuerst (grau gestrichelt)
klassifiziert = set()
for gruppe in (BODEN_KANTEN, ENTHESE_KANTEN, ACHILLES_KANTEN, TALUS_KANTEN):
    for (i, j) in gruppe:
        klassifiziert.add(frozenset((i, j)))
erste = True
for stelle in range(len(RANDWEG)):
    i, j = RANDWEG[stelle], RANDWEG[(stelle + 1) % len(RANDWEG)]
    if frozenset((i, j)) not in klassifiziert:
        kante(i, j, "#8a97a5", 2.4, "--",
              "Rand ohne Lastgruppe" if erste else None)
        erste = False
for (i, j) in BODEN_KANTEN:
    kante(i, j, "#1f77b4", 4.2, "-",
          "Boden" if (i, j) == BODEN_KANTEN[0] else None)
for (i, j) in ENTHESE_KANTEN:            # liegt auf der Bodenfacette (9,10)
    kante(i, j, "#c1440e", 2.0, "-", "Enthese (Faszie)")
for (i, j) in ACHILLES_KANTEN:
    kante(i, j, "#2a7d2a", 4.2, "-", "Achillessehne")
for (i, j) in TALUS_KANTEN:
    kante(i, j, "#7f2fbf", 4.2, "-",
          "Talus-Gelenkfläche" if (i, j) == TALUS_KANTEN[0] else None)

for k in TALUS_LAGERKNOTEN:              # festes Lager als Ring
    ax.plot(kn[k, 0], kn[k, 1], "o", ms=19, mfc="none", mec="#7f2fbf",
            mew=2.4, zorder=5)
for k in range(len(HANDNETZ_KNOTEN)):
    ax.plot(kn[k, 0], kn[k, 1], "o", color="#3d4b59", ms=12, zorder=6)
    ax.text(kn[k, 0], kn[k, 1], str(k), ha="center", va="center",
            fontsize=8, color="white", zorder=7)

ax.plot([44, 54], [-6, -6], "-", color="#333", lw=2)
ax.text(49, -8.5, "10 mm", ha="center", va="top", fontsize=9, color="#333")
ax.text(2, 38.5, "hinten (Ferse)", fontsize=9, color="#555")
ax.text(52, 38.5, "vorn (Zehen)", fontsize=9, color="#555")
ax.legend(loc="lower left", fontsize=7.5, framealpha=0.95)
ax.set_xlim(-14, 72); ax.set_ylim(-19, 42)
ax.set_aspect("equal"); ax.axis("off")
plt.tight_layout()
plt.show()
Abbildung 11.1: Das Handnetz, auf dem dieses Kapitel rechnet: 13 Knoten (nummeriert), 13 Dreiecke, zwei innere Knoten 8 (links) und 11 (rechts) als Fächerzentren. Die Randkanten sind nach der anatomischen Einteilung des Randes gefärbt: blau die drei Bodenkontaktfacetten (0,7), (0,9), (9,10), rot die Faszien-Enthese (9,10) am vorderen unteren Abhang (sie liegt auf der Bodenfacette und ist deshalb doppelt gezeichnet), grün der Achillessehnen-Ansatz (12,7) am hinteren Rand, violett die beiden Talus-Facetten (4,5) und (5,6). Die violetten Ringe um die Knoten 4, 5 und 6 markieren das feste Talus-Lager. Grau gestrichelt die Randkanten ohne Lastgruppe, darunter das Stück (6,12) über dem Achilles-Ansatz. Sagittalschnitt, x nach vorn, y nach oben; Maßstab 10 mm.

Deutung. Drei Gruppen tragen Kraft. Die Bodenreaktion drückt von unten gegen die drei plantaren Facetten (0,7), (0,9) und (9,10); die Plantarfaszie zieht an der Enthese (9,10) am vorderen unteren Abhang — dieselbe Kante trägt also beide Lasten; die Achillessehne zieht am hinteren Rand an der Facette (12,7), deren Endpunkte die Knoten 12 (oben) und 7 (unten) sind. Belastet werden dadurch die Knoten 0, 7, 9, 10 und 12. Festgehalten sind die drei Talus-Knoten 4, 5 und 6 an der oberen Gelenkfläche — dort sitzt die Ferse am Sprungbein. Diese beiden Mengen berühren sich nicht: kein belasteter Knoten ist zugleich gefesselt. Alle übrigen Knoten — die beiden inneren Fächerzentren 8 und 11 sowie die Randknoten 1, 2 und 3 an der Fersenspitze — sind frei und lastfrei; sie bewegen sich allein, weil ihre Nachbarn über die Matrix an ihnen ziehen. Auf dieser Einteilung beruht alles Weitere.

11.3 Die Systemmatrix: dreizehn Elementmatrizen an ihren Plätzen

Jeder Knoten des Handnetzes trägt zwei Freiheitsgrade — eine Verschiebung in x- und eine in y-Richtung. Dreizehn Knoten sind also 26 Freiheitsgrade, und die Systemmatrix ist entsprechend \(26 \times 26\). Die Assemblierung aus Kapitel 10 macht daraus keine Geheimwissenschaft: Sie legt jede der dreizehn \(6 \times 6\)-Elementmatrizen an die Zeilen und Spalten ihrer drei Knoten und addiert überlappende Einträge. Ein einziger Aufruf genügt dafür — mehr steht im folgenden Block nicht:

Die 26×26-Systemmatrix des Handnetzes assemblieren 
HAND_K = assembliere(HANDNETZ_KNOTEN, HANDNETZ_DREIECKE,
                     E_KORTIKAL, NU_KNOCHEN, TIEFE)

Wie groß die entstandene Matrix ist, wie viele ihrer Plätze überhaupt einen Wert tragen und wie der \(2\times 2\)-Block des inneren Knotens 8 aussieht, steht in Tabelle 11.1.

Tabelle 11.1: Größe, Füllgrad und der Diagonalblock des inneren Knotens 8 in der assemblierten Systemmatrix des Handnetzes.
Die assemblierte Systemmatrix K Wert
Größe 26 × 26 Einträge
davon besetzt (ungleich null) 252 von 676
Zeile 16 (Knoten 8, x), Spalten 16 und 17 1530,2 / -194,0 MN/m
Zeile 17 (Knoten 8, y), Spalten 16 und 17 -194,0 / 1475,2 MN/m

Deutung. Von den \(26 \times 26\) Plätzen der Matrix trägt nur ein knappes Drittel einen Wert (Tabelle 11.1) — sie ist eine dünn besetzte Matrix. Der Grund ist geometrisch: Zwei Knoten stehen nur dann gemeinsam in einem Matrixeintrag, wenn sie sich ein Dreieck teilen. Der abgedruckte \(2 \times 2\)-Block gehört zum inneren Knoten 8, der als Zentrum des linken Fächers zu neun Dreiecken gehört; sein Diagonaleintrag ist die Summe der Beiträge aus allen neun und deshalb groß. Ein Randknoten, der nur zu zwei Dreiecken gehört, bekäme einen viel kleineren Diagonalwert. Die Blockstruktur der Matrix ist also ein direktes Abbild der Netztopologie: dicht, wo viele Dreiecke zusammenlaufen (die inneren Knoten 8 und 11), dünn am Rand. Abbildung 11.2 zeigt das ganze Besetzungsmuster als Bild.

Code 
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

K = np.array(HAND_K)
betrag = np.abs(K)
maske = betrag.copy()
maske[maske == 0.0] = np.nan          # exakte Nullen weiss lassen

fig, ax = plt.subplots(figsize=(6.0, 5.4))
vmin = betrag[betrag > 0].min()
vmax = betrag.max()
bild = ax.imshow(maske, cmap="magma_r", norm=LogNorm(vmin=vmin, vmax=vmax))
ax.set_xticks(range(0, 26, 2))
ax.set_yticks(range(0, 26, 2))
ax.set_xticklabels([f"K{k}" for k in range(13)], fontsize=7)
ax.set_yticklabels([f"K{k}" for k in range(13)], fontsize=7)
ax.set_xlabel("Freiheitsgrad (x je Knoten)")
ax.set_ylabel("Freiheitsgrad (x je Knoten)")
for g in range(0, 27, 2):
    ax.axhline(g - 0.5, color="#cccccc", lw=0.3)
    ax.axvline(g - 0.5, color="#cccccc", lw=0.3)
fig.colorbar(bild, ax=ax, fraction=0.046, pad=0.03,
             label="Betrag des Eintrags (N/m, log)")
plt.tight_layout()
plt.show()
Abbildung 11.2: Die 26×26-Systemmatrix des Handnetzes als Besetzungsbild: je dunkler, desto größer der Betrag des Eintrags (logarithmische Farbe). Zwei Freiheitsgrade je Knoten, in Knotenreihenfolge 0…12. Die dichten Blöcke an den Freiheitsgraden 16/17 und 22/23 sind die inneren Knoten 8 und 11 des Doppelfächers, die zu neun bzw. sechs Dreiecken gehören; die dünnen Ränder sind die Randknoten. Weiß = exakt null (die beiden Knoten teilen kein Dreieck). Sagittalschnitt-Netz, Freiheitsgrade x/y je Knoten.

Deutung unter dem Bild. Das Muster ist symmetrisch zur Diagonalen — die Systemmatrix ist symmetrisch, weil die Elementmatrizen es sind (Kapitel 10). Die beiden auffälligen dichten Kreuze liegen bei den Freiheitsgraden der inneren Knoten: Sie „sehen” viele Nachbarn und füllen deshalb ihre Zeilen und Spalten. Diese Assemblierung läuft in der Animation Dreieck für Dreieck ab und geht dort gleich weiter bis zum fertig reduzierten System: Am Knopf füllt sich das Muster sichtbar, Phase für Phase.

11.4 Der Lastvektor: drei Kräfte, auf die Ränder verteilt

Die Matrix beschreibt, wie steif das Netz ist. Was fehlt, ist die rechte Seite des Gleichungssystems — die Last. Für die mittlere Standphase stehen die drei Kräfte im Kräftefahrplan aus Kapitel 2 (Tabelle 2.2): die Bodenreaktion mit \(2{,}5\cdot\)KG nach oben, der Achillessehnenzug mit \(2{,}0\cdot\)KG nach hinten-oben, der Faszienzug mit \(0{,}6\cdot\)KG nach vorn-unten. Keine davon greift an einem einzelnen Punkt an — jede verteilt sich über ihre Randfacetten. Der folgende Block verteilt jede Gesamtkraft längenproportional auf die Endknoten ihrer Randfacetten, je zur Hälfte, und addiert die drei Beiträge zu einem 26er-Lastvektor. Das ist dieselbe konsistente Knotenlast, die Kapitel 9 als natürlichen Randterm hergeleitet hat.

Eine Gesamtkraft auf ihre Randfacetten verteilen und die drei Kräfte zum Lastvektor addieren 
def _facettenlaengen(kanten):
    """Ergaenzt jede Kante (i, j) um ihre Laenge L -> (i, j, L) in Metern."""
    facetten = []
    for (i, j) in kanten:
        (xa, ya) = HANDNETZ_KNOTEN[i]
        (xb, yb) = HANDNETZ_KNOTEN[j]
        facetten.append((i, j, math.hypot(xb - xa, yb - ya)))
    return facetten

BODEN_FAC = _facettenlaengen(BODEN_KANTEN)
ENTHESE_FAC = _facettenlaengen(ENTHESE_KANTEN)
ACHILLES_FAC = _facettenlaengen(ACHILLES_KANTEN)

def verteilte_last(n_knoten, facetten, kraft_xy):
    """Verteilt eine Gesamtkraft (fx, fy) laengenproportional auf die Facetten,
    je zur Haelfte auf ihre Endknoten. Summe der Knotenkraefte = Gesamtkraft."""
    f = [0.0] * (2 * n_knoten)
    laenge_gesamt = sum(fac[2] for fac in facetten)
    fx, fy = kraft_xy
    for (i, j, laenge) in facetten:
        anteil = 0.5 * laenge / laenge_gesamt
        f[2 * i] = f[2 * i] + anteil * fx
        f[2 * i + 1] = f[2 * i + 1] + anteil * fy
        f[2 * j] = f[2 * j] + anteil * fx
        f[2 * j + 1] = f[2 * j + 1] + anteil * fy
    return f

def addiere(a, b):
    return [a[k] + b[k] for k in range(len(a))]

def lastvektor_moment(index):
    """Der 26er-Kanon-Lastvektor eines Moments: die drei Kraefte je Gruppe
    (Faktor*KG*Richtung) verteilt und addiert."""
    _name, f_boden, f_achilles, f_faszie = MOMENTE[index]
    n = len(HANDNETZ_KNOTEN)
    b_kraft = (f_boden * KG * RICHTUNG_BODEN[0], f_boden * KG * RICHTUNG_BODEN[1])
    a_kraft = (f_achilles * KG * RICHTUNG_ACHILLES[0],
               f_achilles * KG * RICHTUNG_ACHILLES[1])
    fa_kraft = (f_faszie * KG * RICHTUNG_FASZIE[0],
                f_faszie * KG * RICHTUNG_FASZIE[1])
    f = verteilte_last(n, BODEN_FAC, b_kraft)
    f = addiere(f, verteilte_last(n, ACHILLES_FAC, a_kraft))
    f = addiere(f, verteilte_last(n, ENTHESE_FAC, fa_kraft))
    return f

f_m2 = lastvektor_moment(2)                 # mittlere Standphase

Tabelle 11.2 zeigt den fertigen Vektor Knoten für Knoten.

Tabelle 11.2: Der Lastvektor der mittleren Standphase (M2): die 26 Knotenkräfte in Newton, je Knoten waagerecht und senkrecht. Positives \(F_x\) zeigt nach vorn, positives \(F_y\) nach oben.
Knoten \(F_x\) (N) \(F_y\) (N)
0 0,0 716,1
1 0,0 0,0
2 0,0 0,0
3 0,0 0,0
4 0,0 0,0
5 0,0 0,0
6 0,0 0,0
7 -220,8 1025,4
8 0,0 0,0
9 217,5 558,3
10 217,5 165,4
11 0,0 0,0
12 -220,8 702,1
größte Knotenkraft auf einem Lagerknoten 0,0

Deutung unter der Tabelle. Nur an fünf Knoten steht eine Kraft: an 0 (reiner Boden), an 7 (wo sich Boden- und Achilles-Facette treffen und ihre Kräfte addieren), an 9 und 10 am vorderen Abhang (wo Boden- und Enthese-Facette zusammenfallen und sich Bodendruck und Faszienzug mischen) sowie an 12, dem oberen Endpunkt der Achilles-Facette, wo allein der Sehnenzug wirkt. Alle anderen Knoten tragen exakt null — sie liegen im Inneren oder auf lastfreien Rändern; die vordere Fersenspitze (Knoten 1 bis 3) ist im Moment M2 unbelastet. Man sieht auch die Richtungen wieder: am reinen Bodenknoten 0 zeigt die Kraft senkrecht nach oben (\(F_y > 0\), \(F_x = 0\)), an den beiden Achilles-Knoten 7 und 12 nach hinten-oben (\(F_x < 0\)), an 9 und 10 zieht der Faszienzug nach vorn (\(F_x > 0\)).

Die letzte Zeile der Tabelle ist die wichtigste für den nächsten Schritt: Auf den drei Talus-Knoten 4, 5 und 6 steht exakt null. Die Achillessehne zieht an der Facette (12, 7) am hinteren Rand, und deren beide Endknoten liegen beide unterhalb der Talus-Gelenkfläche. Last- und Lagerknoten sind damit getrennte Mengen: Jeder Krafteintrag sitzt auf einem Freiheitsgrad, der sich auch bewegen darf. Abbildung 11.3 stellt denselben Vektor als Pfeile über dem Netz dar.

Code 
import numpy as np
import matplotlib.pyplot as plt

kn = np.array(HANDNETZ_KNOTEN) * 1000.0     # m -> mm
f = np.array(lastvektor_moment(2)).reshape(-1, 2)
fig, ax = plt.subplots(figsize=(6.6, 4.4))
for dr in HANDNETZ_DREIECKE:
    ax.fill(kn[list(dr), 0], kn[list(dr), 1], color="#f0e2d0",
            ec="#a9825f", lw=1.2)
skala = 0.010                               # mm je Newton
for k in range(len(HANDNETZ_KNOTEN)):
    fx, fy = f[k]
    if abs(fx) > 1e-9 or abs(fy) > 1e-9:
        if k in (7, 12):
            farbe = "#2a7d2a"               # Achilles (12 rein, 7 mit Boden)
        elif k in (9, 10):
            farbe = "#c1440e"               # Enthese/Faszie (und Boden)
        else:
            farbe = "#1f77b4"               # Boden
        ax.annotate("", xy=(kn[k, 0] + fx * skala, kn[k, 1] + fy * skala),
                    xytext=(kn[k, 0], kn[k, 1]),
                    arrowprops=dict(arrowstyle="->", color=farbe, lw=2.0))
    if k in TALUS_LAGERKNOTEN:              # Lagerknoten: ohne Kraft
        ax.plot(kn[k, 0], kn[k, 1], "o", ms=17, mfc="none", mec="#7f2fbf",
                mew=2.2)
    ax.plot(kn[k, 0], kn[k, 1], "o", color="#5b6b7a", ms=11)
    ax.text(kn[k, 0], kn[k, 1], str(k), ha="center", va="center",
            fontsize=8, color="white")
# Massstab 10 mm
ax.plot([44, 54], [-6, -6], "-", color="#333", lw=2)
ax.text(49, -8.5, "10 mm", ha="center", va="top", fontsize=9, color="#333")
ax.text(2, 37, "hinten (Ferse)", fontsize=9, color="#555")
ax.text(52, 37, "vorn (Zehen)", fontsize=9, color="#555")
ax.set_xlim(-14, 72); ax.set_ylim(-11, 42)
ax.set_aspect("equal"); ax.axis("off")
plt.tight_layout()
plt.show()
Abbildung 11.3: Der Lastvektor der mittleren Standphase (M2) über dem Handnetz: an jedem belasteten Knoten der resultierende Kraftpfeil, Pfeillänge proportional zur Knotenkraft. Blau der reine Bodendruck (Knoten 0) nach oben, grün die Knoten der Achilles-Facette (12 rein, 7 zusammen mit dem Bodendruck) nach hinten-oben, rot die Knoten 9 und 10, an denen sich Bodendruck und Faszienzug mischen. Die drei Talus-Lagerknoten 4, 5 und 6 sind violett umringt und tragen keinen Pfeil — auf ihnen steht exakt null. Die inneren Knoten 8 und 11 und die vordere Fersenspitze bleiben im Moment M2 ebenfalls lastfrei. Sagittalschnitt, x nach vorn, y nach oben; Maßstab 10 mm.

Deutung. Der Lastvektor ist die vollständige rechte Seite des Gleichungssystems. Er enthält keine erfundene Randbedingung — jede der drei Kräfte sitzt an ihrer anatomischen Randgruppe, mit Größe und Richtung aus dem Kräftefahrplan (Tabelle 2.2). Damit ist die Frage „was zieht und drückt an der Ferse?” für diesen Schrittmoment beantwortet; offen bleibt nur noch, was die Ferse festhält.

11.5 Das Lager: die Talus-Knoten festhalten

Ohne Lager hat das Gleichungssystem keine eindeutige Lösung — die Ferse würde unter den Kräften einfach davonfliegen (das war die Kernbotschaft von Kapitel 5). Dieses Modell hält die Ferse an der Talus-Gelenkfläche fest: die Knoten 4, 5 und 6 in beide Richtungen. Der Grund ist nicht, dass ein Gelenk starr wäre — es ist keins. Der Grund ist, dass das Modell die Kraft von oben gar nicht kennt: Das Sprungbein sitzt dort auf und drückt mit einer Kraft, die niemand vorgibt, sondern die sich erst aus dem Gleichgewicht ergibt. Wer den Ort festhält, bekommt diese Kraft als Ergebnis heraus — dazu unten der Abschnitt über die Lagerreaktion. Was die Wahl kostet und welche Bauformen es sonst noch gibt, ordnet Kapitel 12 ein. In Matrix und Lastvektor heißt das, ihre sechs Freiheitsgrade aus dem System auszubauen — die zugehörigen Zeilen und Spalten werden gestrichen, weil ihre Verschiebung von vornherein null ist. Der folgende Block trennt diesen Ausbau (reduziere) vom Lösen (loese_randwert), damit man das kleinere System gleich einzeln anschauen kann.

Die festen Freiheitsgrade ausbauen und das kleinere System lösen 
def reduziere(K, kraefte, feste_liste):
    """Baut die festen Freiheitsgrade aus. Rueckgabe: die Liste der freien
    Freiheitsgrade, die reduzierte Matrix Kr und der reduzierte Lastvektor fr."""
    feste = set(feste_liste)
    frei = []
    for d in range(len(K)):
        if d not in feste:
            frei.append(d)
    Kr = []
    fr = []
    for i in frei:
        neue_zeile = []
        for j in frei:
            neue_zeile.append(K[i][j])
        Kr.append(neue_zeile)
        fr.append(kraefte[i])
    return frei, Kr, fr

def loese_randwert(K, kraefte, feste_liste):
    """Loest K*u = f mit fester Lagerung: die festen Freiheitsgrade werden
    aus Matrix und Lastvektor ausgebaut, das kleinere System geloest, die
    Loesung wieder auf die volle Groesse gesetzt (feste Knoten = 0)."""
    frei, Kr, fr = reduziere(K, kraefte, feste_liste)
    loesung = gauss_loese(Kr, fr)
    u = [0.0] * len(K)
    for stelle in range(len(frei)):
        u[frei[stelle]] = loesung[stelle]
    return u

HAND_FESTE = feste_dofs(TALUS_LAGERKNOTEN)
u_m2 = loese_randwert(HAND_K, lastvektor_moment(2), HAND_FESTE)

Wie viele Freiheitsgrade der Lagereinbau kostet, hält Tabelle 11.3 fest.

Tabelle 11.3: Das Talus-Lager in Zahlen: welche sechs Freiheitsgrade festgehalten werden und wie viele danach noch zu lösen sind.
Der Lagereinbau Wert
feste Freiheitsgrade (Talus-Lager, Knoten 4, 5, 6) 8, 9, 10, 11, 12, 13
Freiheitsgrade gesamt 26
frei nach dem Lagereinbau 20

Deutung unter der Tabelle. Die sechs festen Freiheitsgrade gehören zu den drei Talus-Knoten (je x und y). Nach ihrem Ausbau bleiben von 26 noch 20 freie Freiheitsgrade — das ist das System, das die Gauß-Elimination aus Kapitel 5 tatsächlich löst. Der Lagereinbau ist damit kein Zusatz, sondern erst das, was aus der schwebenden Matrix ein lösbares Problem macht: Die Lagerreaktion an den Talus-Knoten hält die drei aufgebrachten Kräfte im Gleichgewicht. Abbildung 11.4 zeigt, welche Knoten festgehalten werden.

Code 
import numpy as np
import matplotlib.pyplot as plt

kn = np.array(HANDNETZ_KNOTEN) * 1000.0
f_lager = lastvektor_moment(2)
belastet = {k for k in range(len(HANDNETZ_KNOTEN))
            if abs(f_lager[2 * k]) > 1e-9 or abs(f_lager[2 * k + 1]) > 1e-9}
fig, ax = plt.subplots(figsize=(6.6, 4.4))
for dr in HANDNETZ_DREIECKE:
    ax.fill(kn[list(dr), 0], kn[list(dr), 1], color="#f0e2d0",
            ec="#a9825f", lw=1.2)
for k in range(len(HANDNETZ_KNOTEN)):
    if k in TALUS_LAGERKNOTEN:
        farbe = "#c1440e"
    elif k in belastet:
        farbe = "#1f77b4"
    else:
        farbe = "#8a97a5"
    ax.plot(kn[k, 0], kn[k, 1], "o", color=farbe, ms=13)
    ax.text(kn[k, 0], kn[k, 1], str(k), ha="center", va="center",
            fontsize=8, color="white")
for k in TALUS_LAGERKNOTEN:                 # kleines Fesselsymbol
    ax.plot([kn[k, 0] - 2, kn[k, 0] + 2], [kn[k, 1] + 3.2, kn[k, 1] + 3.2],
            "-", color="#c1440e", lw=2)
    for dx in (-1.4, 0.0, 1.4):
        ax.plot([kn[k, 0] + dx, kn[k, 0] + dx - 1.0],
                [kn[k, 1] + 3.2, kn[k, 1] + 4.6], "-", color="#c1440e", lw=1.2)
ax.text(28, 39, "Talus-Lager (fest)", fontsize=9, color="#c1440e")
ax.text(-13.5, 21, "belastet (blau):\n0, 7, 9, 10, 12", ha="left",
        va="center", fontsize=8, color="#1f77b4")
ax.plot([44, 54], [-6, -6], "-", color="#333", lw=2)
ax.text(49, -8.5, "10 mm", ha="center", va="top", fontsize=9, color="#333")
ax.text(2, 37, "hinten (Ferse)", fontsize=9, color="#555")
ax.text(52, 37, "vorn (Zehen)", fontsize=9, color="#555")
ax.set_xlim(-14, 72); ax.set_ylim(-11, 44)
ax.set_aspect("equal"); ax.axis("off")
plt.tight_layout()
plt.show()
Abbildung 11.4: Der Lagereinbau am Handnetz: die drei Talus-Knoten 4, 5, 6 (rot, mit Fesselsymbol) sind in beide Richtungen festgehalten — ihre sechs Freiheitsgrade werden aus dem System ausgebaut. Die belasteten Randknoten 0, 7, 9, 10 und 12 (blau) tragen die drei Kräfte des Moments M2; keiner von ihnen ist zugleich gefesselt. Die freien, lastfreien Knoten (grau) bewegen sich nur über ihre Nachbarn. Aus 26 Freiheitsgraden werden so 20 freie. Sagittalschnitt, x nach vorn, y nach oben; Maßstab 10 mm.

11.5.1 Warum erst die volle Matrix, dann das Streichen

Man könnte fragen, warum die Rechnung erst alle 26 Freiheitsgrade aufbaut und sechs davon anschließend wieder wegwirft, statt sie von vornherein wegzulassen. Drei Gründe sprechen dafür, und keiner davon ist Bequemlichkeit.

Erstens ist die Assemblierung reine Netz- und Materialbuchhaltung. Eine Elementmatrix weiß nur, wie ihre drei Knoten einander Steifigkeit geben — sie kennt Geometrie und Material, sonst nichts. Ob einer dieser Knoten festgehalten wird, ist keine Eigenschaft des Elements, sondern des gestellten Problems. Diese Trennung hält das Verfahren sauber: \(K\) ist die Struktur, \(f\) und die Lager sind der Belastungsfall.

Zweitens braucht man die Zeilen und Spalten der festen Knoten zweimal. Beim Streichen selbst: Weil die Verschiebung dort null ist, fallen ihre Spalten ersatzlos weg — bei einem verschobenen Lager täten sie das nicht, sondern wanderten mit umgekehrtem Vorzeichen als bekannte Anteile auf die rechte Seite. Beides ist dieselbe Bauform von Randbedingung: Der Ort ist vorgegeben, einmal als null, einmal als Wert. Kapitel 12 gibt den Bauformen ihre Namen und stellt ihnen die vorgegebene Kraft und die federnde Bettung zur Seite. Und danach für die Lagerreaktion: Sie ist \(K \cdot u - f\), ausgewertet an genau den gestrichenen Freiheitsgraden. Die gestrichenen Zeilen sind also kein Abfall, sondern der Weg zu der Kraft, mit der das Talus-Gelenk hält.

Drittens zahlt sich die volle Matrix praktisch aus. Dieselbe \(K\) trägt alle fünf Schrittmomente: Solange das Lager dasselbe bleibt, wechselt nur \(f_r\), und assembliert wird trotzdem nur einmal. Sie trägt auch jede andere Lagerwahl — dann ändert sich mit der Liste der freien Freiheitsgrade zwar auch \(K_r\), aber \(K\) selbst bleibt unberührt. Das Streichen ist Buchhaltung, keine neue Rechnung.

11.5.2 Die Randbedingungen in der vollen Matrix

Das Streichen lässt sich in zwei Bewegungen zerlegen, und die erste kann man an der vollen Matrix sehen. Ein festgehaltener Freiheitsgrad sagt zweierlei: Seine Verschiebung ist null, und seine Gleichung wird nicht gelöst, sondern durch die Lagerreaktion erfüllt. Beides zusammen heißt: Die Spalte dieses Freiheitsgrads darf verschwinden (sie würde mit der Verschiebung null multipliziert), und seine Zeile darf verschwinden (sie ist keine Gleichung mehr, die eine Unbekannte bestimmt). Man kann das ausführen, ohne die Matrix kleiner zu machen: Man setzt in jeder betroffenen Zeile und jeder betroffenen Spalte alle Einträge auf null, schreibt eine 1 auf die Diagonale und nullt den zugehörigen Eintrag des Lastvektors. Danach steht dort die triviale Gleichung \(1 \cdot u = 0\) — genau das, was die Lagerbedingung verlangt.

Der folgende Block macht diesen Einbau; wie viele Einträge dabei verschwinden, zählt er gleich mit.

Die Lager in die volle Matrix einbauen, ohne sie zu verkleinern 
def integriere_randbedingungen(K, kraefte, feste_liste):
    """Baut die Lager in die VOLLE Matrix ein, ohne sie zu verkleinern:
    Zeile und Spalte jedes festen Freiheitsgrads werden genullt, auf die
    Diagonale kommt eine 1, der Lasteintrag wird null."""
    K_neu = [zeile[:] for zeile in K]
    f_neu = kraefte[:]
    for d in feste_liste:
        for j in range(len(K)):
            K_neu[d][j] = 0.0
            K_neu[j][d] = 0.0
        K_neu[d][d] = 1.0
        f_neu[d] = 0.0
    return K_neu, f_neu

def zaehle_besetzt(matrix):
    anzahl = 0
    for zeile in matrix:
        for eintrag in zeile:
            if eintrag != 0.0:
                anzahl = anzahl + 1
    return anzahl

K_BC, F_BC = integriere_randbedingungen(HAND_K, lastvektor_moment(2),
                                        HAND_FESTE)

gefallen = 0
for i in range(len(HAND_K)):
    for j in range(len(HAND_K)):
        if HAND_K[i][j] != 0.0 and K_BC[i][j] == 0.0:
            gefallen = gefallen + 1
F_VOR_EINBAU = lastvektor_moment(2)
groesste = 0.0
for d in range(len(HAND_K)):
    unterschied = abs(F_BC[d] - F_VOR_EINBAU[d])
    if unterschied > groesste:
        groesste = unterschied

Die Bilanz des Einbaus steht in Tabelle 11.4.

Tabelle 11.4: Was der Einbau der sechs Lager-Freiheitsgrade in der vollen Matrix bewirkt: verschwundene Kopplungen, neue Diagonal-Einsen und ein unveränderter Lastvektor.
Der Einbau der Randbedingungen Wert
besetzt vor dem Einbau 252
besetzt nach dem Einbau (davon 6 Diagonal-Einsen) 182
auf null gesetzte Koppel-Einträge 70
größte Änderung im Lastvektor 0,0 N

Deutung unter der Tabelle. Ein knappes Drittel der besetzten Einträge fällt weg, und sechs neue Einsen kommen auf der Diagonale hinzu. Verschwunden ist genau das, was die drei festgehaltenen Talus-Knoten mit ihren Nachbarn verband — die Zahl der genullten Kopplungen steht in der Tabelle. Der Lastvektor ändert sich dabei um keinen Newton: Er hatte an den sechs Lager-Freiheitsgraden ohnehin nichts stehen, weil kein belasteter Knoten festgehalten ist. Abbildung 11.5 zeigt die Matrix nach diesem Einbau; rot markiert sind die sechs Zeilen und Spalten der festen Freiheitsgrade.

Code 
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

Kbc = np.abs(np.array(K_BC))
maske_bc = Kbc.copy()
maske_bc[maske_bc == 0.0] = np.nan
fest = sorted(HAND_FESTE)

fig, ax = plt.subplots(figsize=(6.0, 5.4))
# rote Hinterlegung der betroffenen Zeilen und Spalten
for d in fest:
    ax.axhspan(d - 0.5, d + 0.5, color="#f7c9c0", zorder=0)
    ax.axvspan(d - 0.5, d + 0.5, color="#f7c9c0", zorder=0)
vmin = Kbc[Kbc > 0].min()
bild = ax.imshow(maske_bc, cmap="Greys", norm=LogNorm(vmin=vmin, vmax=Kbc.max()),
                 zorder=2)
for d in fest:                       # die neue 1 auf der Diagonalen
    ax.add_patch(plt.Rectangle((d - 0.5, d - 0.5), 1, 1, fill=True,
                               facecolor="#c1440e", edgecolor="#7a2a08",
                               lw=0.6, zorder=3))
ax.set_xticks(range(0, 26, 2)); ax.set_yticks(range(0, 26, 2))
ax.set_xticklabels([f"K{k}" for k in range(13)], fontsize=7)
ax.set_yticklabels([f"K{k}" for k in range(13)], fontsize=7)
ax.set_xlabel("Freiheitsgrad (x je Knoten)")
ax.set_ylabel("Freiheitsgrad (x je Knoten)")
for g in range(0, 27, 2):
    ax.axhline(g - 0.5, color="#cccccc", lw=0.3, zorder=4)
    ax.axvline(g - 0.5, color="#cccccc", lw=0.3, zorder=4)
ax.set_xlim(-0.5, 25.5); ax.set_ylim(25.5, -0.5)
fig.colorbar(bild, ax=ax, fraction=0.046, pad=0.03,
             label="Betrag des Eintrags (N/m, log)")
plt.tight_layout()
plt.show()
Abbildung 11.5: Die volle 26×26-Systemmatrix nach dem Einbau der Randbedingungen: dieselbe Matrix wie Abbildung 11.2, aber Zeile und Spalte jedes der sechs festgehaltenen Freiheitsgrade (Talus-Knoten 4, 5, 6, je x und y) sind genullt und rot hinterlegt; auf ihrer Diagonale steht eine 1 (rotes Kästchen). Graustufe = Betrag des Eintrags (logarithmisch), weiß = exakt null. Man sieht das rote Kreuz aus sechs leeren Zeilen und sechs leeren Spalten — genau diese können im nächsten Schritt wegfallen, ohne dass sich an den übrigen Einträgen etwas ändert. Sagittalschnitt-Netz, Freiheitsgrade x/y je Knoten.

Deutung unter dem Bild. Das rote Kreuz liegt an den Freiheitsgraden 8 bis 13 — das sind die Knoten 4, 5 und 6 mit je x und y. In diesen Zeilen und Spalten steht außer der Diagonal-Eins nichts mehr; sie tragen keine Information über die Verformung des Netzes. Alles außerhalb des Kreuzes ist unverändert: Der Einbau hat keinen einzigen der übrigen Einträge angefasst. Damit ist sichtbar, was das Streichen im nächsten Schritt tut — es entfernt genau die rot markierten Zeilen und Spalten und lässt den Rest, wie er ist.

Der folgende Block führt die zweite Bewegung aus — das eigentliche Streichen, ohne zu lösen — und zeigt damit genau das, was die Gauß-Elimination als Nächstes bekommt: die reduzierte Matrix \(K_r\), den reduzierten Lastvektor \(f_r\) und die Liste der verbliebenen Freiheitsgrade.

Das Streichen ausführen und die verbliebenen Freiheitsgrade benennen 
def dof_name(d):
    """Beschriftung eines Freiheitsgrads: Knotennummer und Richtung."""
    return "K%d %s" % (d // 2, "x" if d % 2 == 0 else "y")

F_VOLL_M2 = lastvektor_moment(2)
FREI_M2, K_R, F_R = reduziere(HAND_K, F_VOLL_M2, HAND_FESTE)
MARKEN_R = []
for d in FREI_M2:
    MARKEN_R.append(dof_name(d))

besetzt_r = 0
for zeile in K_R:
    for eintrag in zeile:
        if eintrag != 0.0:
            besetzt_r = besetzt_r + 1

summe_gestrichen = 0.0
for d in HAND_FESTE:
    summe_gestrichen = summe_gestrichen + abs(F_VOLL_M2[d])

Tabelle 11.5 nennt Größe und Füllgrad des reduzierten Systems und listet die verbliebenen Freiheitsgrade; Tabelle 11.6 zeigt daneben, welche von ihnen überhaupt eine äußere Kraft tragen.

Tabelle 11.5: Das System, das die Gauß-Elimination tatsächlich bekommt: 20 Gleichungen für 20 Unbekannte, und nicht ein einziges Newton geht beim Ausbau verloren.
Das reduzierte System \(K_r \cdot u = f_r\) Wert
Größe von \(K_r\) 20 × 20
Einträge in \(f_r\) 20
in \(K_r\) besetzt (ungleich null) 176 von 400
verbliebene Freiheitsgrade (4, 5, 6 fehlen) K0 x, K0 y, K1 x, K1 y, K2 x, K2 y, K3 x, K3 y, K7 x, K7 y, K8 x, K8 y, K9 x, K9 y, K10 x, K10 y, K11 x, K11 y, K12 x, K12 y
beim Ausbau gestrichene Lasteinträge 0,0 N
Tabelle 11.6: Die besetzten Einträge des reduzierten Lastvektors \(f_r\) — die neun Knotenkräfte aus Tabelle 11.2, die auf freien Freiheitsgraden sitzen.
Freiheitsgrad Kraft (N)
K0 y 716,1
K7 x -220,8
K7 y 1025,4
K9 x 217,5
K9 y 558,3
K10 x 217,5
K10 y 165,4
K12 x -220,8
K12 y 702,1

Deutung unter den Tabellen. Übrig bleiben die zehn Knoten 0, 1, 2, 3, 7, 8, 9, 10, 11 und 12 mit je zwei Richtungen — 20 Gleichungen für 20 Unbekannte. Das ist derselbe Sprung wie in Kapitel 5, nur größer: dort drei Unbekannte einer Federkette, hier zwanzig eines Netzes; das Verfahren dazwischen ist unverändert. Die besetzten Einträge der rechten Seite sind genau die Knotenkräfte aus dem Lastvektor-Abschnitt, keine einzige ist unterwegs verloren gegangen. Die letzte Zeile von Tabelle 11.5 sagt es als Zahl: Beim Ausbau wird kein Newton gestrichen. Die volle Achilleskraft, die volle Bodenreaktion und der volle Faszienzug stehen sämtlich auf freien Freiheitsgraden und wirken damit auch verformend. Abbildung 11.6 zeigt beides zusammen als ein Bild: links das Besetzungsmuster von \(K_r\), rechts die Balken von \(f_r\).

Code 
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

Kr_bild = np.abs(np.array(K_R))
maske_r = Kr_bild.copy()
maske_r[maske_r == 0.0] = np.nan          # exakte Nullen weiss lassen
fr_bild = np.array(F_R)
luecke = 7.5                              # zwischen K3 y und K7 x

fig, (axK, axf) = plt.subplots(
    1, 2, figsize=(7.0, 5.0), layout="constrained",
    gridspec_kw={"width_ratios": [7.0, 2.2]})

bild = axK.imshow(maske_r, cmap="magma_r",
                  norm=LogNorm(vmin=Kr_bild[Kr_bild > 0].min(),
                               vmax=Kr_bild.max()))
axK.set_xticks(range(20)); axK.set_yticks(range(20))
axK.set_xticklabels(MARKEN_R, fontsize=6, rotation=90)
axK.set_yticklabels(MARKEN_R, fontsize=6)
axK.set_xlabel("freier Freiheitsgrad")
axK.set_ylabel("freier Freiheitsgrad")
for g in range(0, 21, 2):
    axK.axhline(g - 0.5, color="#cccccc", lw=0.3)
    axK.axvline(g - 0.5, color="#cccccc", lw=0.3)
axK.axhline(luecke, color="#c1440e", lw=1.4)
axK.axvline(luecke, color="#c1440e", lw=1.4)
axK.set_title("K$_r$ (20 × 20) — rote Linie: hier fehlen\ndie Knoten 4, 5, 6",
              fontsize=9)

axf.barh(range(20), fr_bild,
         color=["#c1440e" if w < 0 else "#1f77b4" for w in fr_bild])
axf.axvline(0.0, color="#555", lw=0.8)
axf.axhline(luecke, color="#c1440e", lw=1.4)
axf.set_ylim(19.5, -0.5)
axf.set_yticks(range(20)); axf.set_yticklabels([])
axf.set_xlabel("Knotenkraft (N)")
axf.tick_params(labelsize=7)
axf.set_title("f$_r$ (20 Einträge)", fontsize=9)

fig.colorbar(bild, ax=[axK, axf], fraction=0.046, pad=0.03,
             label="Betrag des Eintrags in K$_r$ (N/m, log)")
plt.show()
Abbildung 11.6: Das reduzierte System \(K_r \cdot u = f_r\) des Handnetzes für die mittlere Standphase (M2), das die Gauß-Elimination tatsächlich löst. Links die 20×20-Matrix \(K_r\) als Besetzungsbild in derselben Machart wie Abbildung 11.2 (je dunkler, desto größer der Betrag, logarithmische Farbe; weiß = exakt null). Rechts der reduzierte Lastvektor \(f_r\) in Newton als Balken auf denselben 20 Zeilen (blau = positiver Eintrag, also nach vorn beziehungsweise nach oben; rot = negativer, also nach hinten beziehungsweise nach unten). Zeilen und Spalten sind mit Knotennummer und Richtung beschriftet; die rote Linie markiert die Lücke, an der die sechs Freiheitsgrade der festgehaltenen Talus-Knoten 4, 5 und 6 herausgefallen sind. Sagittalschnitt-Netz, Freiheitsgrade x/y je Knoten.

Deutung unter dem Bild. Das linke Muster hat 20 Zeilen und 20 Spalten, und die rote Linie zeigt, wo die Lücke sitzt: zwischen dem letzten Freiheitsgrad des Knotens 3 und dem ersten des Knotens 7 sind die sechs Freiheitsgrade der Knoten 4, 5 und 6 herausgefallen. Weil immer Zeile und Spalte desselben Freiheitsgrads gestrichen werden, bleibt \(K_r\) symmetrisch — das Bild ist wieder spiegelbildlich zur Diagonalen, so wie die volle Matrix es war. Striche man nur die Zeilen, bliebe nicht einmal eine quadratische Matrix übrig — und ohne quadratische Matrix gibt es kein lösbares Gleichungssystem. Der Balken rechts ist die ganze rechte Seite; die Zahlen dazu stehen in Tabelle 11.6. Fünf der zehn verbliebenen Knoten tragen eine Kraft: Knoten 0 trägt reinen Bodendruck nach oben, Knoten 7 Boden und Achilles addiert, die Knoten 9 und 10 den Bodendruck gemischt mit dem Faszienzug nach vorn, Knoten 12 den reinen Sehnenzug nach hinten-oben. Die übrigen elf Einträge sind null. Mehr als diese beiden Objekte bekommt der Löser nicht.

Wer das System nicht als Muster, sondern Zahl für Zahl sehen will, findet es in der folgenden Tabelle: derselbe Zusammenhang, nur ausgeschrieben. Der folgende Block setzt sie — er rechnet nichts, er rundet auf Meganewton je Meter, schreibt exakte Nullen als Punkt und die Diagonale rot und hängt rechts, durch einen Strich abgesetzt, den zugehörigen Eintrag von \(f_r\) in Newton an. Zeile \(i\) der Tabelle ist damit die vollständige \(i\)-te Gleichung des Systems.

Ein Helfer, der das System als scrollbare Zahlentabelle setzt 
def zahlentabelle(matrix, marken, rechte_spalte, rechte_kopf):
    """Scrollbare HTML-Tabelle einer Matrix in MN/m: exakte Nullen als Punkt,
    Diagonale rot, rechts abgesetzt der Lastvektor in Newton."""
    zeilen = ['<div style="overflow:auto; max-height:22em; border:1px solid #ccc; '
              'font-size:11px">', '<table style="border-collapse:collapse">']
    kopf = "<tr><th></th>" + "".join(
        "<th style='padding:2px 4px'>%s</th>" % marken[j]
        for j in range(len(marken)))
    zeilen.append(kopf + "<th style='padding:2px 6px;border-left:2px solid #333'>"
                  + rechte_kopf + "</th></tr>")
    for i in range(len(marken)):
        zellen = ["<th style='padding:2px 4px'>%s</th>" % marken[i]]
        for j in range(len(marken)):
            wert = matrix[i][j] / 1e6      # N/m -> MN/m
            if wert == 0.0:
                zellen.append("<td style='padding:2px 4px;color:#bbb'>·</td>")
            else:
                farbe = "#b00" if i == j else "#333"
                zellen.append("<td style='padding:2px 4px;color:%s'>%s</td>"
                              % (farbe, komma(wert)))
        zellen.append("<td style='padding:2px 6px;border-left:2px solid #333;"
                      "color:#0a6'>%s</td>" % komma(rechte_spalte[i]))
        zeilen.append("<tr>" + "".join(zellen) + "</tr>")
    zeilen.append("</table></div>")
    return "\n".join(zeilen)
K0 x K0 y K1 x K1 y K2 x K2 y K3 x K3 y K7 x K7 y K8 x K8 y K9 x K9 y K10 x K10 y K11 x K11 y K12 x K12 y f [N]
K0 x 446,7 27,3 · · · · · · -154,1 138,7 -26,3 -184,9 -266,2 18,9 · · · · · · 0,0
K0 y 27,3 506,8 · · · · · · 123,3 -262,0 -184,9 -217,8 34,3 -27,0 · · · · · · 716,1
K1 x · · 1190,0 195,4 -239,9 -59,5 · · · · -259,4 124,1 · · -186,2 -119,3 -504,6 -140,7 · · 0,0
K1 y · · 195,4 1534,8 -44,1 -25,9 · · · · 124,1 -18,0 · · -134,7 -449,4 -140,7 -1041,4 · · 0,0
K2 x · · -239,9 -44,1 402,9 52,7 -191,6 -151,2 · · · · · · · · 28,6 142,6 · · 0,0
K2 y · · -59,5 -25,9 52,7 700,2 -135,8 -359,5 · · · · · · · · 142,6 -314,8 · · 0,0
K3 x · · · · -191,6 -135,8 417,4 58,6 · · · · · · · · -144,8 32,2 · · 0,0
K3 y · · · · -151,2 -359,5 58,6 451,8 · · · · · · · · 32,2 256,9 · · 0,0
K7 x -154,1 123,3 · · · · · · 628,3 -34,8 -382,4 -47,2 · · · · · · -91,8 -41,3 -220,8
K7 y 138,7 -262,0 · · · · · · -34,8 634,8 -47,2 4,6 · · · · · · -56,7 -377,4 1025,4
K8 x -26,3 -184,9 -259,4 124,1 · · · · -382,4 -47,2 1530,2 -194,0 53,8 137,2 -378,4 100,8 -152,1 -90,3 -378,1 167,6 0,0
K8 y -184,9 -217,8 124,1 -18,0 · · · · -47,2 4,6 -194,0 1475,2 137,2 -268,4 100,8 -345,4 -90,3 127,5 167,6 -169,5 0,0
K9 x -266,2 34,3 · · · · · · · · 53,8 137,2 628,5 25,8 -416,1 -197,3 · · · · 217,5
K9 y 18,9 -27,0 · · · · · · · · 137,2 -268,4 25,8 592,3 -181,9 -296,9 · · · · 558,3
K10 x · · -186,2 -134,7 · · · · · · -378,4 100,8 -416,1 -181,9 980,6 215,8 · · · · 217,5
K10 y · · -119,3 -449,4 · · · · · · 100,8 -345,4 -197,3 -296,9 215,8 1091,7 · · · · 165,4
K11 x · · -504,6 -140,7 28,6 142,6 -144,8 32,2 · · -152,1 -90,3 · · · · 1620,9 -31,5 · · 0,0
K11 y · · -140,7 -1041,4 142,6 -314,8 32,2 256,9 · · -90,3 127,5 · · · · -31,5 2641,0 · · 0,0
K12 x · · · · · · · · -91,8 -56,7 -378,1 167,6 · · · · · · 707,8 13,5 -220,8
K12 y · · · · · · · · -41,3 -377,4 167,6 -169,5 · · · · · · 13,5 703,9 702,1

Deutung. Zeilen und Spalten tragen Knotennummer und Richtung, und die Knoten 4, 5 und 6 fehlen — die Lücke lässt sich in der Kopfzeile abzählen. Die abgesetzte rechte Spalte ist die Last. Wo sie null ist, sagt die Zeile „an diesem Freiheitsgrad greift keine äußere Kraft an” — die Verschiebung dort ist trotzdem nicht null, weil die Nachbarknoten über die Matrixeinträge der Zeile ziehen und drücken.

11.6 Eine Zahl für Zug, Druck und Schub

Gleich entsteht das erste eingefärbte Feld dieses Buches, das du selbst ausgerechnet hast. Wie man ein Feldbild liest — Farbe, Skala, Orientierung — steht in Kapitel 1; hier geht es um die Größe, die eingefärbt wird.

In jedem Dreieck stehen nach der Lösung drei Spannungen: \(\sigma_{xx}\) in x-Richtung, \(\sigma_{yy}\) in y-Richtung und die Schubspannung \(\tau_{xy}\) (Kapitel 6). Drei Zahlen lassen sich schlecht als eine Farbe malen. Die Von-Mises-Vergleichsspannung fasst sie zu einer zusammen:

\[\sigma_v = \sqrt{\sigma_{xx}^2 - \sigma_{xx}\,\sigma_{yy} + \sigma_{yy}^2 + 3\,\tau_{xy}^2}.\]

Sie ist groß, wenn das Material stark verzerrt wird — gleich ob durch Zug, Druck oder Schub —, und sie ist null, wenn das Material von allen Seiten gleich stark gedrückt wird, wie ein Stein tief unter Wasser. Das macht sie zum bequemen Übersichtsmaß, und deshalb färbt dieses Buch seine Felder damit ein.

Sie hat aber einen blinden Fleck, und er ist für dieses Buch wichtig: Sie unterscheidet Zug nicht von Druck. Ein Dreieck, das mit \(5\) MPa gezogen wird, und eines, das mit \(5\) MPa gedrückt wird, bekommen dieselbe Farbe — setze \(\sigma_{xx} = +5\) oder \(-5\) in die Formel ein, das Quadrat löscht das Vorzeichen aus. Für die Frage, wo Knochen wächst, ist das Vorzeichen aber nicht egal: Ob ein Sporn unter Zug oder unter Druck entsteht, ist genau der Streit, auf den dieses Buch zuläuft. Der folgende Block schreibt die Formel als Funktion auf — mehr steht nicht darin; die Umrechnung auf Kilopascal ist die Einheit, in der das Buch Spannungen nennt.

Die Von-Mises-Vergleichsspannung als Funktion 
def von_mises_kpa(sxx, syy, txy):
    """Von-Mises-Vergleichsspannung (ebener Spannungszustand), in kPa."""
    vm = math.sqrt(sxx * sxx - sxx * syy + syy * syy + 3.0 * txy * txy)
    return vm / 1000.0

Abbildung 11.7 setzt drei Musterzustände in diese Funktion ein.

Code 
import numpy as np
import matplotlib.pyplot as plt

PROBEN = [("Zug\n$\\sigma_{xx} = +5$ MPa", (5e6, 0.0, 0.0), 1),
          ("Druck\n$\\sigma_{xx} = -5$ MPa", (-5e6, 0.0, 0.0), -1),
          ("Schub\n$\\tau_{xy} = 2{,}89$ MPa", (0.0, 0.0, 5e6 / 3 ** 0.5), 0)]

VM_MAX_MPA = 10.0                                    # Skalenende der Abbildung
norm = plt.Normalize(vmin=0.0, vmax=VM_MAX_MPA)
cmap = plt.get_cmap("turbo")

fig, ax = plt.subplots(figsize=(6.8, 3.2))
for stelle, (name, (sxx, syy, txy), pfeil) in enumerate(PROBEN):
    vm = von_mises_kpa(sxx, syy, txy) / 1000.0          # kPa -> MPa
    x0 = stelle * 3.0
    ax.add_patch(plt.Rectangle((x0, 0), 1.8, 1.8, facecolor=cmap(norm(vm)),
                               edgecolor="#222", lw=1.2))
    if pfeil != 0:                       # Zug oder Druck: waagerechtes Paar
        for seite, richtung in ((x0 - 0.15, -1), (x0 + 1.95, +1)):
            ax.annotate("", xy=(seite + pfeil * richtung * 0.55, 0.9),
                        xytext=(seite, 0.9),
                        arrowprops=dict(arrowstyle="->", color="#222", lw=1.8))
    else:                                # Schub: zwei Paare, im Kreis herum
        ax.annotate("", xy=(x0 + 1.5, 2.05), xytext=(x0 + 0.3, 2.05),
                    arrowprops=dict(arrowstyle="->", color="#222", lw=1.8))
        ax.annotate("", xy=(x0 + 0.3, -0.25), xytext=(x0 + 1.5, -0.25),
                    arrowprops=dict(arrowstyle="->", color="#222", lw=1.8))
        ax.annotate("", xy=(x0 + 2.05, 1.5), xytext=(x0 + 2.05, 0.3),
                    arrowprops=dict(arrowstyle="->", color="#222", lw=1.8))
        ax.annotate("", xy=(x0 - 0.25, 0.3), xytext=(x0 - 0.25, 1.5),
                    arrowprops=dict(arrowstyle="->", color="#222", lw=1.8))
    ax.text(x0 + 0.9, 2.55, name, ha="center", va="bottom", fontsize=9)
    ax.text(x0 + 0.9, -0.95, "Von Mises " + komma(vm, 2) + " MPa",
            ha="center", va="top", fontsize=9, color="#333")
ax.set_xlim(-1.2, 8.4); ax.set_ylim(-1.9, 3.6)
ax.set_aspect("equal"); ax.axis("off")
fig.colorbar(plt.cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax,
             fraction=0.046, pad=0.03, label="Von-Mises (MPa)")
plt.tight_layout()
plt.show()
Abbildung 11.7: Der blinde Fleck der Von-Mises-Spannung an drei einachsigen Musterzuständen (schematisch, ohne Ortsbezug): links reiner Zug, in der Mitte reiner Druck, rechts reiner Schub. Die Pfeile zeigen die Belastungsrichtung, die Füllfarbe ist die Von-Mises-Spannung auf der festen Spannungsskala des Buches, hier in Megapascal gelesen (turbo, 0–10 MPa). Zug und Druck sind entgegengesetzt und tragen doch dieselbe Farbe; der Schub ist so gewählt, dass er dieselbe Vergleichsspannung ergibt. Die Zahl unter jedem Feld ist das Ergebnis der Formel.

Deutung. Alle drei Felder haben dieselbe Farbe, obwohl das linke Material auseinandergezogen, das mittlere zusammengedrückt und das rechte nur verschert wird. Die Von-Mises-Farbe sagt also, wo viel passiert, nicht was passiert. Bis Kapitel 15 die Vorzeichen zurückholt — dort werden Zug und Druck an jeder Stelle getrennt —, ist sie das Übersichtsmaß dieses Buches, und mehr nicht.

11.7 Die Lösung: das erste Feldbild aus reiner Handrechnung

Jetzt ist alles beisammen: Matrix, Lastvektor, Lager. Der letzte Schritt übersetzt die gelöste Verschiebung in eine Spannung. Aus der konstanten Dehnung jedes P1-Dreiecks (Kapitel 10) folgt über das Material seine konstante Spannung, und daraus die eben eingeführte Vergleichsspannung. Der folgende Block rechnet für M2 die Verschiebung der Bodenknoten und die dreizehn Von-Mises-Werte, fasst alles in einer Funktion rechne_moment zusammen, die später auch die anderen Momente bedient, und hält am Ende die höchste Spannung jedes der fünf Momente fest.

Von der gelösten Verschiebung zur Von-Mises-Spannung je Dreieck 
def mittlere_uy_mm(u, knoten_liste):
    """Mittlere senkrechte Verschiebung einer Knotengruppe (mm, plus = oben)."""
    summe = sum(u[2 * k + 1] for k in knoten_liste)
    return summe / len(knoten_liste) * 1000.0

def max_verschiebung_mm(u, n_knoten):
    gross = 0.0
    for k in range(n_knoten):
        betrag = math.hypot(u[2 * k], u[2 * k + 1])
        gross = max(gross, betrag)
    return gross * 1000.0

def element_spannung(ecken, u, dr, e_modul, nu):
    """Konstante Spannung [sxx, syy, txy] (Pa) eines P1-Dreiecks."""
    b, c, _ = zeltsteigungen(ecken)
    ux = [u[2 * dr[0]], u[2 * dr[1]], u[2 * dr[2]]]
    uy = [u[2 * dr[0] + 1], u[2 * dr[1] + 1], u[2 * dr[2] + 1]]
    exx = b[0] * ux[0] + b[1] * ux[1] + b[2] * ux[2]
    eyy = c[0] * uy[0] + c[1] * uy[1] + c[2] * uy[2]
    gam = (c[0] * ux[0] + c[1] * ux[1] + c[2] * ux[2]
           + b[0] * uy[0] + b[1] * uy[1] + b[2] * uy[2])
    fk = e_modul / (1.0 - nu * nu)
    sxx = fk * (exx + nu * eyy)
    syy = fk * (nu * exx + eyy)
    txy = e_modul / (2.0 * (1.0 + nu)) * gam
    return sxx, syy, txy

def rechne_moment(index):
    """Die ganze Handrechnung eines Moments: Lastvektor, Lagereinbau, Loesung,
    Von-Mises je Element. Rueckgabe als dict."""
    f = lastvektor_moment(index)
    u = loese_randwert(HAND_K, f, HAND_FESTE)
    vm = []
    for dr in HANDNETZ_DREIECKE:
        ecken = [HANDNETZ_KNOTEN[dr[0]], HANDNETZ_KNOTEN[dr[1]],
                 HANDNETZ_KNOTEN[dr[2]]]
        sxx, syy, txy = element_spannung(ecken, u, dr, E_KORTIKAL, NU_KNOCHEN)
        vm.append(von_mises_kpa(sxx, syy, txy))
    return {"u": u, "von_mises_kpa": vm,
            "boden_senkung_mm": mittlere_uy_mm(u, [0, 7, 9, 10]),
            "max_u_mm": max_verschiebung_mm(u, len(HANDNETZ_KNOTEN))}

erg = rechne_moment(2)
SPITZEN_KPA = [max(rechne_moment(k)["von_mises_kpa"])
               for k in range(len(MOMENTE))]

Tabelle 11.7 fasst die Verschiebungen zusammen, Tabelle 11.8 nennt die dreizehn Von-Mises-Werte Dreieck für Dreieck.

Tabelle 11.7: Die Verformung des Handnetzes in der mittleren Standphase — Ausgaben dieses Modells, keine festen Zahlen des Buches.
Die Lösung des Moments M2 Wert
mittlere senkrechte Verschiebung der Bodenknoten (+ = nach oben) 0,008022 mm
größter Verschiebungsbetrag 0,0158 mm
Tabelle 11.8: Die Von-Mises-Vergleichsspannung je Dreieck des Handnetzes im Moment M2. Jedes P1-Dreieck trägt genau einen konstanten Wert.
Dreieck Knoten Von-Mises (kPa)
0 (7, 0, 8) 3290,8
1 (0, 9, 8) 4452,2
2 (9, 10, 8) 3138,9
3 (10, 1, 8) 1734,5
4 (5, 6, 8) 4211,5
5 (6, 12, 8) 8840,0
6 (12, 7, 8) 5627,8
7 (1, 2, 11) 1379,2
8 (2, 3, 11) 566,0
9 (3, 4, 11) 1032,9
10 (4, 5, 11) 1131,6
11 (1, 11, 8) 800,5
12 (5, 8, 11) 1923,4

Deutung unter den Tabellen. Die mittlere senkrechte Verschiebung der Bodenknoten ist mit 0,0080 mm winzig und zeigt nach oben: Die Bodenreaktion drückt die plantare Kante gegen das feste Talus-Lager darüber — es ist eine Zusammendrückung, kein Absinken. Kortikaler Knochen ist steif, und das Handnetz ist grob. Wichtiger sind die dreizehn Von-Mises-Werte aus Tabelle 11.8: Der kleinste steht im lastfernen vorderen Fächer (Dreieck 2, 3, 11), der größte im obersten Dreieck am hinteren Rand (Knoten 6, 12, 8), und zwischen beiden liegt mehr als ein Faktor zehn. Die Spitze sitzt also dort, wo die Achillessehne zieht — genauer: in dem Dreieck zwischen dem festgehaltenen Talus-Knoten 6 und dem gezogenen Achilles-Knoten 12, das den kurzen Weg vom Krafteintrag zum Lager überbrückt. Nicht, weil das Bild es so malt, sondern weil die Rechnung es so ausgibt.

WarnungDas grobe Handnetz liefert nur eine Karikatur der Spannung

Naheliegende Vermutung nach dem ersten bunten Feldbild: „Jetzt haben wir die Spannung in der Ferse.”

Warum sie naheliegt: Das Bild sieht aus wie ein Messergebnis — dreizehn Farben, eine klare Spitze, eine Zahl in Kilopascal.

Was stattdessen gilt: Dreizehn konstante Elementwerte sind kein belastbares Spannungsfeld. Jedes Dreieck trägt genau einen Wert; die Spitze an der Krafteinleitung ist grob unteraufgelöst, und der Zahlenwert hängt an der Netzfeinheit — der Regler am Kapitelende rechnet genau dieses Problem auf den feineren Netzen und findet dort ein Vielfaches dieser Spitze. Das Handnetz zeigt den Weg — assemblieren, lasten, lagern, lösen —, nicht die belastbare Zahl. Der Weg zur belastbaren Zahl ist die Verfeinerung, und ihre saubere Prüfung folgt in Kapitel 13.

Wie grob, sieht man am Ort der Spitze: Sie sitzt in allen fünf Momenten im selben Dreieck (6, 12, 8) — dem einen Dreieck, das zwischen dem Achilles-Krafteintrag und der Talus-Lagerkante liegt. Ein Netz mit dreizehn Dreiecken hat für den Weg vom Krafteintrag zum Lager eben nur dieses eine Dreieck; wo genau die Spannung ihr Maximum hat und wie hoch es ist, kann es gar nicht auflösen. Der Zahlenwert dort ist eine Größenordnung, keine Messung.

Dass diese Handrechnung mit dem schnellen Browser-Löser übereinstimmt, ist kein Zufall, sondern eine algebraische Identität: Beide diskretisieren dieselbe schwache Form mit denselben P1-Dreiecken. Das Buch prüft bei jedem Bau, dass der sichtbare Lehrcode dieses Kapitels und der schnelle Löser hinter den Reglern dasselbe herausbekommen wie eine dritte, unabhängig geschriebene Rechnung: für alle fünf Momente besser als ein Milliardstel (\(10^{-10}\)), in Verschiebung und Spannung. Wo sie auseinanderliefen, wäre das ein Fehler, kein Spielraum. Abbildung 11.8 zeigt das verformte Netz und das Von-Mises-Feld der mittleren Standphase.

Code 
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.tri import Triangulation

kn = np.array(HANDNETZ_KNOTEN) * 1000.0
tri = Triangulation(kn[:, 0], kn[:, 1], np.array(HANDNETZ_DREIECKE))
res = rechne_moment(2)
vm = np.array(res["von_mises_kpa"])
u = np.array(res["u"]).reshape(-1, 2) * 1000.0     # mm
kn_v = kn + 400.0 * u                               # ueberhoehtes Netz

fig, ax = plt.subplots(figsize=(6.8, 4.8))
bild = ax.tripcolor(tri, facecolors=np.clip(vm, 0, 11000), cmap="turbo",
                    vmin=0, vmax=11000, edgecolors="none")
ax.triplot(kn[:, 0], kn[:, 1], np.array(HANDNETZ_DREIECKE),
           color="#d8c6ac", lw=0.6)
ax.triplot(kn_v[:, 0], kn_v[:, 1], np.array(HANDNETZ_DREIECKE),
           color="#111111", lw=1.1)
cb = fig.colorbar(bild, ax=ax, fraction=0.046, pad=0.03)
cb.set_label("Von-Mises (kPa)")
ax.text(2, 37, "hinten (Ferse)", fontsize=9, color="#555")
ax.text(52, 37, "vorn (Zehen)", fontsize=9, color="#555")
ax.plot([44, 54], [-8, -8], "-", color="#333", lw=2)
ax.text(49, -10.5, "10 mm", ha="center", va="top", fontsize=9, color="#333")
ax.text(2, -10.5, "Verformung 400-fach überhöht", fontsize=8, color="#333")
ax.set_xlim(-14, 72); ax.set_ylim(-13, 42)
ax.set_aspect("equal"); ax.axis("off")
plt.tight_layout()
plt.show()
Abbildung 11.8: Das erste Feldbild aus reiner Handrechnung (Moment M2, mittlere Standphase): Von-Mises je Element (turbo), darüber das 400-fach überhöhte verformte Netz (dunkle Kanten) über dem Ausgangsnetz (hell). Die Skala ist ausdrücklich 0–11 000 kPa benannt — sie ist über alle fünf Momente dieselbe und reicht bis über die höchste von ihnen (die Fersenablösung); die feste Spannungsskala des Buches 0–1000 kPa wäre gesättigt. Werte-Legende rechts, Ortszeile unten. Sagittalschnitt, x nach vorn, y nach oben; Maßstab 10 mm.

11.7.1 Was das Gelenk hält: die Lagerreaktion

Kapitel 2 hat drei Kräfte benannt und eine vierte offen gelassen: die Kraft, mit der das Sprungbein von oben auf die Ferse drückt. Sie wurde nicht vorgegeben, weil niemand sie kennt; sie ergibt sich erst aus dem Gleichgewicht. Genau diese Kraft steht jetzt in den Zeilen, die beim Streichen weggefallen sind.

Der Weg dahin ist eine Zeile Rechnung. Für einen freien Freiheitsgrad gilt nach der Lösung \(K \cdot u = f\) — innere Kräfte gleich äußere. Für einen festgehaltenen gilt das nicht: Dort ist \(K \cdot u\) im Allgemeinen nicht gleich \(f\), und was fehlt, muss das Lager beisteuern. Die Lagerreaktion ist deshalb die Differenz \(K \cdot u - f\), ausgewertet an den festen Freiheitsgraden. Der folgende Block rechnet sie aus.

Die Lagerreaktion an den festgehaltenen Freiheitsgraden 
def lagerreaktion(K, u, kraefte, feste_liste):
    """Die Kraft, die das Lager beisteuert: (K*u - f) an den festen
    Freiheitsgraden. Rueckgabe: dict Freiheitsgrad -> Kraft in Newton."""
    reaktion = {}
    for d in feste_liste:
        summe = 0.0
        for j in range(len(K)):
            summe = summe + K[d][j] * u[j]
        reaktion[d] = summe - kraefte[d]
    return reaktion

R_M2 = lagerreaktion(HAND_K, u_m2, lastvektor_moment(2), HAND_FESTE)

Tabelle 11.9 stellt die drei Lagerknoten und ihre Summe der Summe aller aufgebrachten Kräfte gegenüber.

Tabelle 11.9: Die Lagerreaktion an der Talus-Gelenkfläche in der mittleren Standphase, Knoten für Knoten, und ihre Bilanz gegen die drei aufgebrachten Kräfte.
\(F_x\) (N) \(F_y\) (N)
Reaktion am Knoten 4 277,8 513,4
Reaktion am Knoten 5 -145,1 -642,1
Reaktion am Knoten 6 -126,1 -3038,7
Summe der Lagerreaktionen 6,6 -3167,4
Summe der aufgebrachten Kräfte -6,6 3167,4
Summe aus beidem 0,0 0,0
Betrag der Lagerreaktion 3167,5 N 4,30 · KG

Deutung unter der Tabelle. Die letzte Bilanzzeile ist null — in beiden Richtungen, bis auf Rundungsstellen. Das ist kein Zufall und keine gute Näherung, sondern das Gleichgewicht der ganzen Ferse: Was Boden, Sehne und Faszie an Kraft einbringen, trägt das Talus-Gelenk ab, und mehr Wege hinaus gibt es in diesem Modell nicht. Der Betrag ist die gesuchte vierte Kraft — gut vier Körpergewichte, nach unten in die Ferse gerichtet, obwohl der Boden nur zweieinhalb davon hochdrückt. Der Rest kommt aus der Achillessehne, die von hinten-oben zieht und deren Gegenzug ebenfalls das Gelenk aufnehmen muss.

Auf die einzelnen Knoten verteilt sich diese Kraft dagegen sehr ungleich: Am hinteren Lagerknoten 6, gleich neben dem Sehnenansatz, drückt das Gelenk fast so stark wie in der Summe, am vorderen Knoten 4 zieht es sogar nach oben. Auch das ist Mechanik und kein Fehler: Die Ferse kippt unter dem Sehnenzug gegen das Gelenk, und ein Kippen braucht an einem Ende Druck und am anderen Zug. Wie sich diese Verteilung ändert, wenn das Netz feiner wird und das Gelenk aus vielen Knoten statt aus dreien besteht, ist eine Frage an Kapitel 12; die Summe ändert sich dabei nicht, denn sie hängt nur am Gleichgewicht — der Regler am Kapitelende zeigt es.

11.8 Der ganze Weg als Film

Die vier Schritte des Fahrplans lassen sich am besten laufend zeigen — an dem Gegenstand, an dem sie alle stattfinden: an der Matrix selbst. Die Animation läuft in vier Phasen. Phase 1 baut die \(26 \times 26\)-Matrix Dreieck für Dreieck ein; bei jedem Schritt leuchten die sechs Freiheitsgrade des eben hinzugefügten Dreiecks auf, und das Besetzungsmuster wächst, bis alle dreizehn Dreiecke stehen. Phase 2 füllt den Lastvektor, der rechts neben der Matrix als schmale Spalte steht: Knoten für Knoten, wobei nur die fünf belasteten Knoten überhaupt einen Balken bekommen. Phase 3 baut die Randbedingungen ein — die sechs Zeilen und Spalten der Talus-Freiheitsgrade färben sich rot und nullen sich, eine nach der anderen. Phase 4 streicht sie weg: Die roten Zeilen und Spalten fallen heraus, und die Matrix zieht sich auf die \(20 \times 20\) zusammen, die der Löser bekommt.

11.9 Die fünf Schrittmomente am Regler

Bisher stand ein einziger Moment im Mittelpunkt — die mittlere Standphase. Der Schritt hat aber fünf Momentaufnahmen, jede mit anderer Mischung aus Boden-, Achilles- und Faszienkraft — sie stehen im Kräftefahrplan (Tabelle 2.2). Dieselbe rechne_moment-Funktion bedient alle fünf: Der Umschalter zeigt für jeden Moment das verformte Netz und das Von-Mises-Feld.

WichtigVorhersage-Punkt

Bevor du umschaltest: Bei welchem der fünf Schrittmomente ist die höchste Spannung im Fersenbein am größten — beim Fersenauftritt (viel Boden, wenig Sehne), in der mittleren Standphase oder beim Abstoß (wenig Boden, viel Sehne und Faszie)? Lege dich fest.

Deutung. Die höchste Spitze tritt weder beim Fersenauftritt (M0, 3005 kPa) noch beim Abstoß (M4, 6201 kPa) auf, sondern bei der Fersenablösung (M3) mit 10379 kPa; danach folgen die mittlere Standphase (M2, 8840 kPa) und der Belastungsaufbau (M1, 7005 kPa). Der Ort ist in allen fünf Momenten derselbe: das Dreieck (6, 12, 8) am hinteren Rand. Die Reihenfolge folgt keiner der beiden einfachen Regeln. Nach dem Sehnenzug allein müsste M4 (\(2{,}8\cdot\)KG) über M2 (\(2{,}0\cdot\)KG) liegen — er liegt darunter. Nach dem Bodendruck allein müsste M2 (\(2{,}5\cdot\)KG) über M3 (\(1{,}2\cdot\)KG) liegen — auch das stimmt nicht. Was die Spitze treibt, ist die Überlagerung: M3 verbindet den mit Abstand größten Sehnenzug mit einem noch spürbaren Bodendruck, M4 hat zwar viel Sehne, aber fast keinen Boden mehr, der dagegenhält. Wer allein auf „viel Boden gleich viel Spannung” oder allein auf den Sehnenzug getippt hat, unterschätzt das Zusammenspiel — die Rechnung zeigt beides.

11.10 Netzfeinheit: dasselbe Problem, mehr Dreiecke

Das Handnetz war grob mit Absicht — man sollte jeden Schritt sehen. Der Weg selbst — assemblieren, lasten, lagern, lösen — hängt aber nicht an dreizehn Dreiecken. Der letzte Regler geht ihn noch einmal, und zwar mit genau dem Problem dieses Kapitels: dieselben drei Kräfte der mittleren Standphase auf ihren Randfacetten, dieselbe festgehaltene Talus-Gelenkfläche. Verändert wird allein das Netz — vom Handnetz mit dreizehn Dreiecken zu den drei Buchnetzen aus Kapitel 10, deren feinstes über viertausend hat.

Zwei Zahlen meldet er zu jeder Stufe: die höchste Von-Mises-Spannung und die Summe der Lagerreaktion aus dem vorigen Abschnitt.

WichtigVorhersage-Punkt

Bevor du schiebst: Zwei Fragen, lege dich bei beiden fest. Nähert sich die höchste Spannung mit feinerem Netz einem festen Wert, oder wächst sie weiter? Und: Ändert sich die Kraft, mit der das Talus-Gelenk hält?

Deutung. Die beiden Zahlen verhalten sich vollkommen verschieden, und das ist die Lehre dieses Reglers.

Die Lagerreaktion steht auf allen vier Stufen bis auf die letzte angezeigte Stelle gleich. Sie muss es: Sie ist nichts als das Gleichgewicht — was hinein geht, muss heraus —, und das Gleichgewicht kennt keine Netzfeinheit. Wer wissen will, mit welcher Kraft das Sprungbein auf die Ferse drückt, bekommt sie schon vom Handnetz mit dreizehn Dreiecken richtig.

Die höchste Spannung dagegen wächst mit jeder Verfeinerung weiter, vom Handnetz bis zum feinsten Buchnetz um fast eine Größenordnung. Sie pendelt sich nicht ein. Der Grund ist im Bild zu sehen: Die Spitze rückt mit jeder Stufe näher an das hintere Ende der festgehaltenen Talus-Fläche — an die Stelle also, wo ein fest eingespannter Rand auf einen freien trifft. Dort hat schon das gerechnete Modell keine endliche Spannung; das Netz malt nur, wie weit es dem Unendlichen bisher gefolgt ist. Das Handnetz hatte für diese Stelle ein einziges Dreieck und musste sie deshalb unterschätzen.

Daraus folgt die Regel, mit der Kapitel 13 weiterarbeitet: Eine Größe, die aus einer Bilanz kommt, ist früh belastbar; eine Größe, die an einer Ecke ihr Maximum hat, ist es womöglich nie. Welche Größen dieses Buch deshalb auswertet und wie man das prüft, ist die Frage des nächsten Teils.

Damit schließt sich der Bogen dieses Kapitels. Vier Schritte standen am Anfang, vier Ergebnisse stehen am Ende: Die Systemmatrix war die Netztopologie in Zahlen, dicht an den inneren Knoten und leer zwischen Knoten ohne gemeinsames Dreieck. Der Lastvektor war keine Erfindung, sondern die drei Kräfte des Kräftefahrplans, verteilt auf ihre Randfacetten. Das Lager war kein Zusatz, sondern erst das, was die Aufgabe lösbar machte — und es gab die vierte Kraft zurück, die Kapitel 2 offen gelassen hatte. Die Lösung war ein Feldbild, das den Ort der Spitze richtig andeutet und ihre Höhe nicht.

Übungen

Ü 11.1 (Verstehen). Der Diagonaleintrag der Systemmatrix ist am inneren Knoten 8 groß und an einem Randknoten klein. Begründe das, ohne zu rechnen — zähle, zu wie vielen Dreiecken jeder der beiden Knoten gehört.

Ein Knoten bekommt seinen Diagonaleintrag aus jedem Dreieck, das ihn enthält: Beim Zusammenbauen addiert die Assemblierung die Beiträge aller angrenzenden Elemente auf denselben Platz. Die Zahl der Nachbardreiecke ist damit unmittelbar der Grund für die Größe des Eintrags.

Im Handnetz zählt man am inneren Knoten 8 — dem Zentrum des linken Fächers — neun Dreiecke. Ein typischer Randknoten liegt am Ende einer Kette und gehört nur zu zwei. Neun Beiträge gegen zwei: Der Diagonaleintrag am Knoten 8 muss ein Vielfaches des Randwertes sein, ohne dass man einen einzigen davon ausrechnet.

Die Auszählung über alle dreizehn Knoten zeigt dasselbe Muster: Nur die zwei inneren Knoten 8 und 11 haben hohe Grade (neun und sechs), zwei Randknoten kommen auf drei, alle übrigen neun auf zwei. Deshalb ist das Besetzungsmuster in Abbildung 11.2 dort dicht, wo viele Dreiecke zusammenlaufen, und dünn am Rand — die Matrix ist ein Abbild der Netztopologie.

Ü 11.2 (Verändern). Schalte den Momente-Umschalter durch alle fünf Momente und protokolliere die höchste Von-Mises-Spannung. Sage vorher, ob die Spitze beim Fersenauftritt oder beim Abstoß größer ist, und prüfe deine Vorhersage.

Der Umschalter liefert für M0 bis M4 die höchsten Von-Mises-Werte 3005, 7005, 8840, 10379 und 6201 kPa.

Zur gestellten Frage: Die Spitze beim Abstoß ist gut doppelt so groß wie beim Fersenauftritt. Wer auf den Fersenauftritt getippt hat, weil dort der Aufprall sitzt, wird widerlegt — beim Abstoß zieht die Achillessehne, und dieser Zug ist die größere Belastung.

Wichtiger ist aber, was die Frage nicht fragt: Das Maximum liegt bei keinem der beiden Momente, sondern dazwischen, bei der Fersenablösung. Die Reihe steigt von M0 bis M3 an und fällt zu M4 wieder ab. Wer nur die zwei genannten Momente vergleicht, übersieht den größten Wert des ganzen Schrittes.

Ein Vorbehalt gehört dazu: Das sind Spitzenwerte einzelner Elemente auf einem Netz aus dreizehn Dreiecken. Kapitel 13 zeigt, dass ein solcher Einzelelementwert am Netz hängt und keine belastbare Zahl ist. Die Reihenfolge der fünf Momente ist die Aussage, die hier trägt — nicht die Stellen hinter der führenden Ziffer.

Ü 11.3 (Übertragen). Die Handrechnung und die scikit-fem-Referenz stimmen auf besser als \(10^{-10}\) überein. Begründe, warum das eine algebraische Identität sein muss und nicht bloß eine gute Näherung — welche Zutat teilen beide Löser?

Die geteilte Zutat ist die Diskretisierung: dieselbe schwache Form, dieselben P1-Dreiecke. Damit ist alles festgelegt, was in die Rechnung eingeht — die Ansatzfunktionen, die Elementmatrix, der Lastvektor, der Einbau des Lagers.

Entscheidend ist, dass beim linearen Dreieck nichts genähert werden muss: Die Dehnung ist im Element konstant, also lässt sich das Integral der Elementsteifigkeit in geschlossener Form hinschreiben. Es gibt keine Quadraturregel, bei der zwei Umsetzungen unterschiedlich genau sein könnten. Beide Löser bauen deshalb Zeichen für Zeichen dieselbe Matrix \(K\) und denselben Vektor \(f\) und lösen dasselbe endliche Gleichungssystem \(K \cdot u = f\).

Was übrig bleibt, ist nicht Genauigkeit, sondern Reihenfolge: Die beiden Programme addieren und eliminieren in anderer Folge, und Gleitkommazahlen sind nicht exakt assoziativ. Übrig bleibt Rundungsrauschen — und das ist der Grund, warum die Schranke bei \(10^{-10}\) liegt und nicht bei null.

Daraus folgt auch, was die Übereinstimmung nicht zeigt. Sie sagt nichts darüber, wie gut das Modell die Wirklichkeit trifft: Beide Löser haben denselben Diskretisierungsfehler gegenüber dem Kontinuumsproblem, und der fällt beim Vergleich heraus, statt sichtbar zu werden. Zwei Programme, die dasselbe falsche Netz rechnen, stimmen ebenso perfekt überein. Ob das Netz fein genug ist, ist eine andere Frage — die von Kapitel 13.

Roter Faden

Zurück: Die Elementmatrix und die Assemblierung aus Kapitel 10 werden hier zur vollen Systemmatrix; die Einspannung aus Kapitel 5 wird das Talus-Lager, die natürlichen Randterme aus Kapitel 9 werden die verteilte Last; die Von-Mises-Spannung fasst die drei Spannungen des ebenen Zustands aus Kapitel 6 zu einer Farbe zusammen. Vor: Die Bauformen der Randbedingung ordnet Kapitel 12 systematisch am feinen Netz — das Handnetz hat festes Lager und verteilte Last hier erstmals konkret angewandt, und die Lagerreaktion dieses Kapitels ist die Gelenkkraft, die Kapitel 2 offen gelassen hatte. Dass die Von-Mises-Farbe Zug und Druck nicht trennt, holt Kapitel 15 nach. Wie fein das Netz sein muss und warum ein einzelnes buntes Feldbild noch keine belastbare Zahl ist, prüft Kapitel 13.

Was dieses Kapitel NICHT tut

Es liefert keine belastbare Spannungszahl für die Ferse — dreizehn konstante Elementwerte sind eine Karikatur (siehe Widerlegungskasten). Es prüft nicht die Netzkonvergenz (das ist Kapitel 13) und führt keine neue Randbedingungstheorie ein (Kapitel 12). Es rechnet keine neue feste Zahl: Die Verschiebungs- und Von-Mises-Werte sind Ausgaben dieses Modells; das Material stammt unverändert aus dem Zahlenanhang (Anhang A), die Aufteilung der drei Kräfte aus dem Kräftefahrplan (Tabelle 2.2). Es vernetzt nicht im Browser (gmsh bleibt offline) und mittelt keinen Zeitverlauf — die fünf Momente sind Momentaufnahmen.