10 Dreiecke über die Ferse

In Kapitel 9 wurde die Gleichgewichtsbedingung abgeschwächt: Statt „gleich null in jedem Punkt” verlangt die schwache Form nur noch „gleich null im gewichteten Mittel gegen jede zulässige Testfunktion” — und sie braucht dafür nur erste Ableitungen. Damit passt sie zu geknickten, stückweise geraden Ansätzen. Dieses Kapitel macht diese Ansätze konkret: Die Fersenkontur wird mit P1-Dreieck-Elementen überzogen, jedes Dreieck bekommt seine kleine Elementsteifigkeitsmatrix, und der Zusammenbau (Assemblierung) fügt sie zur großen Matrix \(K\) aus Kapitel 4 und Kapitel 5. Damit schließt sich die Kette „Körper → Netz → Gleichungssystem” — und zwar so, dass man sie nachrechnen kann: auf einem winzigen Netz, das im Buch abgedruckt steht.

Lernziele

Nach diesem Kapitel kannst du

erklären, warum man die Fersenkontur in Dreiecke zerlegt und was ein Knoten, eine Kante und ein Element ist,
an einem einzelnen P1-Dreieck die lineare Ansatzfunktion und ihre konstante Dehnung nachvollziehen,
verstehen, wie die kleine Elementmatrix aus der schwachen Form, dem Material und der Dreiecksgeometrie entsteht und wie der Zusammenbau daraus die große Matrix K macht,
am winzigen Handnetz die Assemblierung zweier Dreiecke Eintrag für Eintrag selbst nachrechnen und sehen, wie sich an ihrer geteilten Kante die Beiträge beider Dreiecke addieren.

10.1 Der Aufhänger: gerade Dreiecke für einen gebogenen Spannungsverlauf

Ein Rechner kennt keine Kurven, nur Zahlen an Punkten. Um die Ferse zu rechnen, wird ihr Gebiet in lauter kleine, gerade begrenzte Dreiecke zerlegt — ein Netz; dieses Zerlegen heißt Vernetzung. Jedes Dreieck rechnet stur linear: eine schräge Ebene, eine konstante Dehnung. Im Knochen ändert sich die Belastung aber weich und gebogen von Ort zu Ort. Wie kann ein Haufen gerader Dreiecke, von denen jedes nur eine einzige Dehnung kennt, einen solchen gebogenen Verlauf nachbilden — und wie viele braucht man dafür?

Bevor das erste Netz zu sehen ist, wird das Rechenwerk dieses Kapitels eingerichtet. In ihm steht alles, was hier gerechnet wird: die Steigungen der Zelt-Ansatzfunktionen, die Materialtabelle des Hookeschen Gesetzes, die Elementmatrix, der Zusammenbau, der Gauß-Löser aus Kapitel 5, dazu das winzige Handnetz und der Dreiecksform-Demonstrator der späteren Abschnitte. Reines Python, ohne eine einzige Rechenbibliothek — jede Zeile ließe sich mit Papier und Bleistift nachvollziehen. Die Zelle ist eingeklappt: Wer wissen will, woher ein Name in den Blöcken weiter unten kommt, klappt sie auf. Material und Scheibendicke holt sie aus dem Zahlenanhang (Anhang A), sie erfindet keine Zahl.

Das Rechenwerk dieses Kapitels — zum Nachlesen aufklappen

# Reines Python: außer der Werteliste des Zahlenanhangs keine Bibliothek,
# keine Abkürzung. Numerisch gleich gehalten mit js/kap10_dreiecke.js
# (window.Kap10Dreiecke.mathe); beide werden gegeneinander und gegen eine
# unabhängige scikit-fem-Rechnung geprüft (Abweichung < 10⁻¹⁰).
import math
import sys

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

def komma(x, stellen=2):
    """Deutsche Kommadarstellung mit Endnullen-Trimmung; Minuszeichen −."""
    s = f"{x:.{stellen}f}"
    if "." in s:
        s = s.rstrip("0").rstrip(".")
    return s.replace("-", "−").replace(".", ",")

def komma_sig(x, n=3):
    """Deutsche Darstellung mit n signifikanten Stellen (ohne e-Notation);
    echte Rundungsreste (< 10⁻¹²) werden als 0 gezeigt, sonst NIE nullgeschnappt."""
    if abs(x) < 1e-12:
        return "0"
    s = f"{x:.{n}g}"
    if "e" in s or "E" in s:
        s = f"{x:.12f}".rstrip("0").rstrip(".")
    return s.replace("-", "−").replace(".", ",")

# --- Ein einzelnes P1-Dreieck ----------------------------------------------
def zeltsteigungen(ecken):
    """Die konstanten Steigungen der drei Zelt-Ansatzfunktionen eines Dreiecks.
    ecken = [(x,y), (x,y), (x,y)]. Rückgabe (b, c, flaeche) mit b[i] = dN_i/dx,
    c[i] = dN_i/dy (je 1/m) und flaeche = Dreiecksfläche (m², positiv). Weil jede
    Ansatzfunktion über dem Dreieck eine schräge Ebene ist, ist ihre Steigung
    KONSTANT — der einzige Geometrie-Baustein, den wir brauchen."""
    (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(f"Dreieck entartet (Fläche null): {ecken}")
    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 3×3-Materialtabelle C des ebenen Spannungszustands (Hooke, Kapitel 6).
    Verknüpft die drei Dehnungen [eps_xx, eps_yy, gamma_xy] mit den drei
    Spannungen [sigma_xx, sigma_yy, tau_xy]. nu koppelt die beiden Richtungen."""
    if e_modul <= 0.0:
        raise ValueError(f"E-Modul muss positiv sein: {e_modul}")
    if not (0.0 <= nu < 0.5):
        raise ValueError(f"Querkontraktion ν muss in [0; 0,5) liegen: {nu}")
    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):
    """Die 3×6-Tabelle B: aus sechs Eckverschiebungen [u0,v0,u1,v1,u2,v2] die drei
    Dehnungen. Zeile 0: eps_xx = Summe b_i·u_i; Zeile 1: eps_yy = Summe c_i·v_i;
    Zeile 2: gamma_xy = Summe c_i·u_i + b_i·v_i."""
    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 6×6-Elementsteifigkeitsmatrix eines P1-Dreiecks: B-transponiert · C · B,
    mal Tiefe mal Fläche. Weil B konstant ist, ist das Integral nur „Wert mal
    Fläche". Eintrag (i, j) ist die Kraft am Freiheitsgrad i, wenn man den
    Freiheitsgrad j um eins verschiebt."""
    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

# --- Zusammenbau (Assemblierung) und Löser ---------------------------------
def leeres_system(dof):
    return [[0.0] * dof for _ in range(dof)]

def assembliere(knoten, dreiecke, e_modul, nu, tiefe):
    """Setzt die große Systemmatrix K aus den kleinen Elementmatrizen zusammen.
    Für jedes Dreieck den 6×6-Steckbrief bauen und jeden seiner 36 Einträge an die
    Freiheitsgrade seiner drei Ecken addieren. Knoten i besetzt Zeile/Spalte 2·i
    (waagerecht) und 2·i+1 (senkrecht). Wo Dreiecke sich einen Knoten teilen,
    summieren sich ihre Beiträge."""
    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 b in range(6):
                K[dofs[a]][dofs[b]] = K[dofs[a]][dofs[b]] + Ke[a][b]
    return K

def gauss_loese(A, rhs):
    """Gauß-Elimination mit Spaltenpivot (wie Kapitel 5), für kleine Systeme."""
    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 lösbar (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 element_dehnung(ecken, eck_u):
    """Die drei Dehnungen eines Dreiecks aus seinen sechs Eckverschiebungen (B·d).
    Ein einziger Wert je Dreieck — die konstante Dehnung des P1-Elements."""
    B = b_tabelle(ecken)
    d = [0.0, 0.0, 0.0]
    for z in range(3):
        summe = 0.0
        for k in range(6):
            summe = summe + B[z][k] * eck_u[k]
        d[z] = summe
    return d

def element_spannung(ecken, eck_u, e_modul, nu):
    """Die drei Spannungen eines Dreiecks: sigma = C · (B · Eckverschiebungen)."""
    d = element_dehnung(ecken, eck_u)
    C = materialtabelle(e_modul, nu)
    s = [0.0, 0.0, 0.0]
    for z in range(3):
        summe = 0.0
        for k in range(3):
            summe = summe + C[z][k] * d[k]
        s[z] = summe
    return s

def von_mises(spannung):
    """Von-Mises-Vergleichsspannung aus [sigma_xx, sigma_yy, tau_xy]."""
    sxx, syy, txy = spannung
    return math.sqrt(sxx * sxx - sxx * syy + syy * syy + 3.0 * txy * txy)

# --- Das Handnetz (13 Knoten, 13 Dreiecke) — winzig, abdruckbar, von Hand --
# Identisch mit dem Handnetz der Referenzrechnung und mit js/kap10_dreiecke.js.
# Doppelfächer: zwei innere Knoten 8 (links) und 11 (rechts), zwischen ihnen
# ein Verbindungsband (Dreiecke (1,11,8) und (5,8,11)). Der plantare Fächerteil
# 0 → 9 → 10 → 1 bildet den vorderen unteren Abhang des Tuber nach: Knoten 9
# (21,8/0,0 mm) ist der tiefste Punkt (medialer Fortsatz), Knoten 10 (25,9/7,1 mm)
# das obere Ende des Abhangs; die Facette 9 → 10 ist die Enthese am realen
# Sporn-Ort. Am hinteren Rand sitzt zwischen den Knoten 6 (12,0/32,0 mm)
# und 7 (1,0/12,0 mm) der Randknoten 12 (4,0/22,0 mm); die beiden Fächerzellen
# (6,12,8) und (12,7,8) füllen den hinteren Rand.
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)]
E_KORTIKAL = kanon.wert_von("E_KORTIKAL")   # Pa, kortikaler Knochen
NU_KNOCHEN = kanon.wert_von("NU")           # Querkontraktion Knochen
TIEFE = kanon.wert_von("DICKE")             # m, Ersatzdicke der Scheibe

# --- Dreiecksform-Demonstrator (schematisch, E = 1, ν = 0,3, ohne Einheit) --
FORM_RAHMEN = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)]  # A,B,C,D fest
FORM_DREIECKE = [(0, 1, 4), (1, 2, 4), (2, 3, 4), (3, 0, 4)]     # um P (Index 4)
FORM_E, FORM_NU, FORM_TIEFE = 1.0, 0.30, 1.0

def form_knoten(spitzigkeit):
    if not (0.0 <= spitzigkeit <= 1.0):
        raise ValueError(f"Spitzigkeit muss in [0; 1] liegen: {spitzigkeit}")
    py = 0.5 - 0.48 * spitzigkeit
    return FORM_RAHMEN + [(0.5, py)]

def dreieck_guete(ecken):
    """Formgüte eines Dreiecks: 1 für gleichseitig, gegen 0 für einen Splitter
    (4·√3·Fläche geteilt durch die Summe der Kantenquadrate)."""
    (xa, ya), (xb, yb), (xc, yc) = ecken
    _, _, flaeche = zeltsteigungen(ecken)
    l2 = ((xb - xa) ** 2 + (yb - ya) ** 2
          + (xc - xb) ** 2 + (yc - yb) ** 2
          + (xa - xc) ** 2 + (ya - yc) ** 2)
    return 4.0 * math.sqrt(3.0) * flaeche / l2

def form_min_guete(spitzigkeit):
    kn = form_knoten(spitzigkeit)
    return min(dreieck_guete([kn[d[0]], kn[d[1]], kn[d[2]]]) for d in FORM_DREIECKE)

def form_max_eintrag(spitzigkeit):
    """Größter Betrag in der Elementmatrix des spitzesten Dreiecks — wächst
    ∝ 1/Fläche, wenn das Dreieck zum Splitter wird."""
    kn = form_knoten(spitzigkeit)
    beste, groesster = None, 0.0
    for d in FORM_DREIECKE:
        ecken = [kn[d[0]], kn[d[1]], kn[d[2]]]
        g = dreieck_guete(ecken)
        Ke = elementmatrix(ecken, FORM_E, FORM_NU, FORM_TIEFE)
        m = max(abs(v) for zeile in Ke for v in zeile)
        if beste is None or g < beste:
            beste, groesster = g, m
    return groesster

def form_reduziert(spitzigkeit):
    """Die reduzierte 2×2-Systemmatrix des einen freien Knotens P (Index 4)."""
    kn = form_knoten(spitzigkeit)
    K = assembliere(kn, FORM_DREIECKE, FORM_E, FORM_NU, FORM_TIEFE)
    return [[K[8][8], K[8][9]], [K[9][8], K[9][9]]]

def form_kondition(spitzigkeit):
    """Konditionszahl der reduzierten 2×2-Matrix (größter/kleinster Eigenwert)."""
    (a, b), (_, d) = form_reduziert(spitzigkeit)
    mitte = (a + d) / 2.0
    radius = math.sqrt(((a - d) / 2.0) ** 2 + b * b)
    lmax, lmin = mitte + radius, mitte - radius
    if lmin <= 0.0:
        raise ValueError("reduzierte Matrix nicht positiv definit")
    return lmax / lmin

def form_senkung(spitzigkeit):
    """Senkung des freien Knotens unter einer Einheitslast nach unten."""
    u = gauss_loese(form_reduziert(spitzigkeit), [0.0, -1.0])
    return math.sqrt(u[0] * u[0] + u[1] * u[1])

Der zweite Teil zeichnet nur. Er lädt die echte 48-Punkt-Kontur des Fersenbeins und die drei Buchnetze, mit denen die Regler dieses Buches rechnen, setzt unter jedes Bild die Orientierungszeile mit Maßstab und holt die aufwendigeren Darstellungen — die Farbkarte der Elementmatrix, das Besetzungsmuster — aus programme/kap10/werkbank_bilder.py, damit das Matplotlib-Beiwerk den Lehrgedanken der Blöcke weiter unten nicht überwuchert.

Die Zeichenwerkzeuge der Bilder dieses Kapitels

# --- Nur Zeichenhilfen (numpy/matplotlib): echte Kontur, echte Buchnetze,
# Orientierung, Maßstab. KEINE Lehr-Mathematik — die steht im Rechenwerk oben.
import json
import pathlib
import sys

import numpy as np
import matplotlib.pyplot as plt

# Zeichenmechanik der sichtbaren Bilder: der Lehrgedanke bleibt in den
# sichtbaren Zellen, das Matplotlib-Beiwerk liegt in programme/kap10
# (Heatmap der Elementmatrix, Besetzungsmuster).
sys.path.insert(0, str(pathlib.Path("programme/kap10").resolve()))
from werkbank_bilder import zeichne_elementmatrix, besetzung_gross

# Echte Kalkaneus-Kontur (48 Punkte, 66 × 42,06 mm) und die drei echten
# Buchnetze (aus derselben Kontur vernetzt).
_kontur = json.loads(
    pathlib.Path("../geometrie/kalkaneus_kontur.json").read_text())
KONTUR = np.array(_kontur["punkte"]) * 1000.0        # m -> mm

BUCHNETZE = {}
for _stufe in ("grob", "mittel", "fein"):
    _d = json.loads(pathlib.Path(f"netze/netz_{_stufe}.json").read_text())
    BUCHNETZE[_stufe] = {
        "knoten": np.array(_d["knoten"]) * 1000.0,   # m -> mm
        "dreiecke": np.array(_d["dreiecke"]),
        "n_knoten": _d["n_knoten"], "n_dreiecke": _d["n_dreiecke"],
    }

def orientierung(ax, mit_massstab=True):
    """„oben ↑"-Pfeil, hinten/vorn und 10-mm-Maßstabsbalken."""
    ax.text(-13, 43, "hinten (Ferse)", fontsize=8, color="#666")
    ax.text(44, 43, "vorn (Zehen)", fontsize=8, color="#666")
    ax.annotate("oben", xy=(-11, 40), xytext=(-11, 26), fontsize=8, color="#444",
                ha="center", arrowprops=dict(arrowstyle="->", color="#444", lw=1.4))
    if mit_massstab:
        ax.plot([54, 64], [-2, -2], color="#333", lw=2.4, solid_capstyle="butt")
        ax.plot([54, 54], [-3.2, -0.8], color="#333", lw=1.2)
        ax.plot([64, 64], [-3.2, -0.8], color="#333", lw=1.2)
        ax.text(59, -5.4, "10 mm", ha="center", fontsize=8, color="#333")

Abbildung 10.1 zeigt dieselbe Fersenkontur einmal grob und einmal fein vernetzt.

Code

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

fig, achsen = plt.subplots(3, 1, figsize=(6.4, 9.6))
for ax, stufe in zip(achsen, ("grob", "mittel", "fein")):
    netz = BUCHNETZE[stufe]
    kn = netz["knoten"]
    tri = Triangulation(kn[:, 0], kn[:, 1], netz["dreiecke"])
    ax.triplot(tri, color="#a9825f", lw=0.5)
    ax.fill(KONTUR[:, 0], KONTUR[:, 1], color="#f0e2d0", ec="#8a6a4a",
            lw=1.4, zorder=0)
    ax.set_title(f"{stufe}: {netz['n_knoten']} Knoten, "
                 f"{netz['n_dreiecke']} Dreiecke", fontsize=10, pad=10)
    orientierung(ax)
    ax.set_xlim(-16, 70); ax.set_ylim(-8, 52)
    ax.set_aspect("equal"); ax.axis("off")
plt.tight_layout()
plt.show()

Abbildung 10.1: Die drei Buchnetze über derselben echten Fersenkontur, untereinander: grob (330 Knoten, 585 Dreiecke), mittel (768 / 1416) und fein (2199 / 4186). Jede Ecke eines Dreiecks ist ein Knoten, jede Dreieckseite eine Kante. Von oben nach unten wächst nur die Zahl der Dreiecke; die feste Polygonkontur des Randes bleibt in allen drei Stufen unverändert — dieselbe Form, feiner aufgelöst. Sagittalschnitt, x nach vorn, y nach oben; Maßstab 10 mm.

Deutung. Es ist immer dieselbe Kontur — nur die Zahl der Dreiecke wächst. Schon das grobe Netz zeichnet die Form erkennbar nach. Der Rand ist dabei in allen drei Stufen dieselbe feste Polygonkontur; die Verfeinerung unterteilt genau deren gerade Segmente in mehr und kleinere Dreiecke, statt eine Krümmung feiner nachzuzeichnen. Was feiner wird, ist also nicht der Umriss, sondern die Dichte der Dreiecke im Inneren — dieselbe Form, feiner aufgelöst. Die Antwort auf die Aufhängerfrage lautet: Jedes Dreieck bleibt stur linear und kennt nur seine eine Dehnung, aber viele nebeneinander setzen daraus einen Verlauf zusammen, der von Dreieck zu Dreieck weiterspringt — je feiner das Netz, desto kleiner die Sprünge und desto glatter der gebogene Verlauf. Wie viele Dreiecke es genau braucht, ist eine eigene Frage — sie gehört zu Kapitel 13.

Ein Netz ist keine Zeichnung

Vermutung: „Feiner gezeichnete Dreiecke sehen glatter aus, also ist das Ergebnis schöner und richtiger.”

Warum sie naheliegt: Das feine Netz in Abbildung 10.1 wirkt sauberer, und im Alltag heißt „höher aufgelöst” meist „besser”.

Was stattdessen gilt: Das Netz ist keine Illustration, sondern die Struktur der Rechnung. Ein Netz, das an der Lasteinleitung zu grob ist, liefert dort falsche Spannungen — auch wenn die Kontur ordentlich aussieht. Und umgekehrt kann ein Netz im ruhigen Inneren grob bleiben, ohne zu schaden. Wo verfeinert werden muss, entscheidet die Physik (wie schnell sich die Spannung örtlich ändert), nicht die Optik. Genau das ist der Grund für die Konvergenzstudie in Kapitel 13 — und der Grund, warum die berechnete Randverschiebung weiter unten am Regler nicht sprunghaft, aber auch nicht völlig ruhig ist.

10.2 Das einzelne Dreieck: drei Knoten, eine schräge Ebene

Der Baustein ist das P1-Dreieck-Element: drei Ecken (Knoten), und dazwischen wird das Verschiebungsfeld linear interpoliert. Zu jedem Knoten gehört eine Ansatzfunktion \(N_i\) — eine schräge Ebene, die im eigenen Knoten den Wert eins hat und in den beiden anderen null. Abbildung 10.2 zeigt so eine Zeltfunktion über einem Dreieck.

Code

import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure(figsize=(8.4, 4.2))
ecken = np.array([[0.0, 0.0], [1.0, 0.0], [0.35, 0.9]])

axl = fig.add_subplot(1, 2, 1)
axl.fill(ecken[:, 0], ecken[:, 1], color="#f0e2d0", ec="#a9825f", lw=1.8)
for i, (x, y) in enumerate(ecken):
    axl.plot(x, y, "o", color="#5b6b7a", ms=9)
    axl.text(x, y - 0.08, f"Knoten {i}", ha="center", va="top", fontsize=9)
axl.text(0.45, 0.3, "lineare\nInterpolation", ha="center", fontsize=9, color="#8a6a4a")
axl.set_xlim(-0.2, 1.2); axl.set_ylim(-0.25, 1.1)
axl.set_aspect("equal"); axl.axis("off")
axl.set_title("drei Knoten", fontsize=10)

axr = fig.add_subplot(1, 2, 2, projection="3d")
X = np.array([ecken[0, 0], ecken[1, 0], ecken[2, 0]])
Y = np.array([ecken[0, 1], ecken[1, 1], ecken[2, 1]])
axr.plot_trisurf(X, Y, np.zeros(3), color="#f0e2d0", alpha=0.8,
                 edgecolor="#a9825f")
axr.plot_trisurf(X, Y, np.array([0.0, 1.0, 0.0]), color="#b3541e",
                 alpha=0.55, edgecolor="#7a2b0e")
axr.plot([ecken[1, 0], ecken[1, 0]], [ecken[1, 1], ecken[1, 1]], [0.0, 1.0],
         "--", color="#7a2b0e", lw=1.2)
axr.plot([ecken[1, 0]], [ecken[1, 1]], [1.0], "o", color="#b3541e", ms=7)
axr.text(ecken[1, 0], ecken[1, 1], 1.12, "N₁ = 1", fontsize=9, color="#7a2b0e")
axr.view_init(elev=22, azim=-135)
axr.set_title("Ansatzfunktion N₁ als schräge Ebene", fontsize=10)
axr.set_zlim(0, 1.2); axr.set_xticks([]); axr.set_yticks([]); axr.set_zticks([])
axr.set_xlabel("x"); axr.set_ylabel("y")
plt.tight_layout()
plt.show()

Abbildung 10.2: Ein P1-Dreieck (schematisch, ohne Einheit). Links: drei Knoten, dazwischen wird linear interpoliert. Rechts: die Ansatzfunktion N₁ des unteren rechten Knotens als schräge Ebene über der grau gezeichneten Dreiecks-Grundfläche — eins im eigenen Knoten (gestrichelte Steiglinie), null in den anderen. Weil die Ebene schräg, aber eben ist, ist ihre Steigung — und damit die Dehnung du/dx im Element — überall gleich: eine konstante Dehnung je Dreieck. Rein schematisch, kein Maßstab.

Deutung. Weil die Ansatzfunktion eine ebene schräge Fläche ist, ist ihre Steigung in \(x\)- und in \(y\)-Richtung überall dieselbe. Setzt man das Verschiebungsfeld aus solchen Ebenen zusammen, ist auch seine Ableitung — die Dehnung — in jedem Dreieck konstant. Diese konstante Steigung ist die einzige Geometriegröße, die wir brauchen; im Rechenwerk heißt sie zeltsteigungen. Tabelle 10.1 führt sie an einem konkreten Dreieck vor: rechtwinklig, mit 20 mm langen Katheten. Dazu wird seine Ecke 1 um 0,01 mm nach vorn geschoben, während die beiden anderen Ecken stehen bleiben — und element_dehnung sagt, welche Dehnung daraus folgt.

Tabelle 10.1: Ein rechtwinkliges P1-Dreieck mit 20 mm langen Katheten: die drei Zeltsteigungen, die Fläche und die Dehnung, die aus einer Verschiebung der Ecke 1 um 0,01 mm in x folgt. Die drei Dehnungswerte gelten für das GANZE Dreieck.

Größe	Ecke 0	Ecke 1	Ecke 2
Steigung \(b_i = \partial N_i/\partial x\) (1/m)	−50	50	0
Steigung \(c_i = \partial N_i/\partial y\) (1/m)	−50	0	50
Fläche des Dreiecks (mm²)	200
Dehnung \(\varepsilon_{xx}\)	0,0005
Dehnung \(\varepsilon_{yy}\)	0
Schub \(\gamma_{xy}\)	0

Deutung unter der Tabelle. Die drei Steigungen sind reine Geometrie (aus den Eckkoordinaten), und die Dehnung, die daraus für eine Eckverschiebung folgt, ist ein Zahlentripel für das ganze Dreieck — nicht ein Feld, das sich innen ändert. Ein feineres Netz braucht man also gerade deshalb, weil ein einzelnes Element die Dehnung nur als Konstante darstellen kann.

Die drei Dehnungen entstehen direkt aus den sechs Eckverschiebungen — jede Ecke \(i\) hat eine waagerechte Verschiebung \(u_i\) und eine senkrechte \(v_i\), gewichtet mit den Steigungen \(b_i\) (in \(x\)) und \(c_i\) (in \(y\)):

\[ \varepsilon_{xx} = \sum_i b_i\,u_i, \qquad \varepsilon_{yy} = \sum_i c_i\,v_i, \qquad \gamma_{xy} = \sum_i \big(c_i\,u_i + b_i\,v_i\big). \]

Das sind genau die drei Verformungen, die Kapitel 6 am Würfelchen eingeführt hat: die Dehnung \(\varepsilon_{xx}\) in \(x\), die Dehnung \(\varepsilon_{yy}\) in \(y\) und die Gleitung \(\gamma_{xy}\), der Winkel, um den der rechte Winkel kippt. Neu ist hier nur, woher sie kommen — nämlich aus den sechs Eckverschiebungen, gewichtet mit reiner Geometrie. Mehr steckt nicht darin: Diese drei Summen sind der ganze Schritt von den Eckverschiebungen zur Dehnung, und weil die Steigungen \(b_i\) und \(c_i\) im Dreieck konstant sind, ist es auch das Ergebnis.

10.3 Das Handnetz: ein winziges Netz zum Nachrechnen von Hand

Damit Elementmatrix und Zusammenbau gleich an einem greifbaren Beispiel entstehen, braucht es zuerst ein konkretes, kleines Netz — klein genug, um es abzudrucken und jeden Schritt von Hand nachzurechnen. Das ist das Handnetz in Abbildung 10.3: dreizehn nummerierte Knoten, dreizehn Dreiecke, eine grobe Karikatur der Fersenform. Genau dieses Netz begleitet das ganze Buch — die Kapitel 11 bis 13 und Kapitel 17 rechnen daran weiter, und die Nummern der Knoten bleiben dabei dieselben.

Zwei seiner Dreiecke sind hervorgehoben, weil die nächsten Abschnitte an ihnen arbeiten: das Dreieck (0, 9, 8) und seine Nachbarin (9, 10, 8). Beide teilen sich die Kante 8 → 9. Die Kante 9 → 10 der zweiten ist die Enthese — die Ansatzfläche der Plantarfaszie, an der Kapitel 11 später zieht und um die es am Ende des Buches geht. Die Knotenkoordinaten sind Millimeter im Sagittalschnitt, gemessen vom Nullpunkt aus Kapitel 1 — hinten unten am Knochen, erste Koordinate nach vorn, zweite nach oben. Deshalb sitzt Knoten 9, der tiefste Punkt des Knochens, genau bei \(y = 0\).

Code

import numpy as np
import matplotlib.pyplot as plt

kn = np.array(HANDNETZ_KNOTEN) * 1000.0
zwei = [(0, 9, 8), (9, 10, 8)]
fig, ax = plt.subplots(figsize=(6.6, 4.6))
for dr in HANDNETZ_DREIECKE:
    ax.fill(kn[list(dr), 0], kn[list(dr), 1], color="#f0e2d0", ec="#a9825f", lw=1.4)
for dr, farbe in zip(zwei, ("#e7c9a0", "#d9a066")):
    ax.fill(kn[list(dr), 0], kn[list(dr), 1], color=farbe, ec="#a9825f", lw=1.4)
ax.plot(kn[[8, 9], 0], kn[[8, 9], 1], color="#7a2b8a", lw=3.4,
        solid_capstyle="round")
for i, (x, y) in enumerate(kn):
    ax.plot(x, y, "o", color="#5b6b7a", ms=13)
    ax.text(x, y, str(i), ha="center", va="center", fontsize=9, color="white")
ax.text(11.5, 4.5, "(0, 9, 8)", fontsize=8, color="#6a4a1a")
ax.text(24.0, 11.5, "(9, 10, 8)", fontsize=8, color="#6a4a1a")
ax.annotate("geteilte Kante 8 → 9",
            xy=((kn[8, 0] + kn[9, 0]) / 2.0, (kn[8, 1] + kn[9, 1]) / 2.0),
            xytext=(0, 27), fontsize=9, color="#7a2b8a",
            arrowprops=dict(arrowstyle="->", color="#7a2b8a", lw=1.2))
orientierung(ax)
ax.set_xlim(-16, 70); ax.set_ylim(-8, 44)
ax.set_aspect("equal"); ax.axis("off")
plt.tight_layout()
plt.show()

Abbildung 10.3: Das Handnetz: dreizehn nummerierte Knoten, dreizehn Dreiecke, eine grobe Karikatur der Fersenform. Hervorgehoben die beiden Dreiecke, an denen unten die Elementmatrix und die Assemblierung entstehen — (0, 9, 8) und (9, 10, 8); violett ihre geteilte Kante 8 → 9. Die Kante 9 → 10 des zweiten Dreiecks ist die Enthese, die Ansatzfläche der Plantarfaszie. Knoten 9 ist der tiefste Punkt des Knochens und liegt deshalb bei y = 0. Sagittalschnitt, x nach vorn, y nach oben; Maßstab 10 mm.

Deutung. Das Handnetz ist bewusst grob — es soll nicht die echte Ferse ersetzen, sondern den Weg zeigen, den die nächsten Abschnitte gehen: erst der Steckbrief eines einzelnen seiner Dreiecke (0, 9, 8), dann der Zusammenbau zweier benachbarter Steckbriefe an ihrer geteilten Kante 8 → 9 und schließlich die Verallgemeinerung auf alle dreizehn. Es ist eine Bühne, an der sich Elementmatrix und Assemblierung an echten Koordinaten zeigen lassen — noch ohne Last und Lager; die vollständige Rechnung dieses Netzes mit den Lasten des Zahlenanhangs und dem Talus-Lager folgt in Kapitel 11.

10.4 Die Elementsteifigkeitsmatrix: Material trifft Geometrie

Die konstante Dehnung wird über das Hookesche Gesetz (Kapitel 6) zur Spannung, und aus Dehnung, Spannung und Fläche entsteht die Elementsteifigkeitsmatrix — der „Steckbrief” des Dreiecks. Wie viele Zeilen und Spalten er hat, sagt die Zahl der Freiheitsgrade: In Kapitel 4 trug jeder Knoten der Federkette einen einzigen, hier trägt jeder Knoten zwei — einen waagerechten und einen senkrechten. Ein Dreieck hat drei Ecken, also sechs Freiheitsgrade, und sein Steckbrief ist eine \(6\times 6\)-Tabelle. Jeden Eintrag kann man ohne jede Formel lesen, in derselben Spaltenlesart wie damals: Eintrag \((i,j)\) ist die Kraft, die am Freiheitsgrad \(i\) entsteht, wenn man allein den Freiheitsgrad \(j\) um eins verschiebt und alle anderen festhält. Jede Spalte ist also die Antwort des Dreiecks — sechs Knotenkräfte — auf eine einzelne Einheits-Eckverschiebung.

Von der schwachen Form zum Steckbrief. Woher kommt dieser Steckbrief? Aus der schwachen Form von Kapitel 9. Dort stand für den Stab \(\int_0^L E\,A\,u'\,v'\,\mathrm{d}x\): die Dehnung der Lösung mal Material mal die Dehnung der Testfunktion, über die Länge aufsummiert. In der Ebene ist es dasselbe, nur dass „Dehnung” jetzt drei Zahlen sind — \(\varepsilon_{xx}\), \(\varepsilon_{yy}\) und die Gleitung \(\gamma_{xy}\) — und „Material” die 3×3-Tabelle \(C\) aus Kapitel 6. Aus dem Produkt zweier Zahlen wird also die Summe über drei Anteile, \(\int \varepsilon^{\top} C\,\varepsilon\,\mathrm{d}A\), und aus dem Stabstück ein Dreieck. Jetzt kommt der Schritt, der alles einfach macht: Setzt man für Lösung und Testfunktion die Zeltfunktionen dieses Dreiecks ein — beide aus demselben Vorrat, das ist das Galerkin-Verfahren aus Kapitel 9 —, dann sind beide Dehnungen im Dreieck konstant. Ein Integral über einen konstanten Integranden aber ist nichts weiter als Wert mal Fläche (mal Dicke, weil die Scheibe eine Tiefe hat). Genau das ist die Formel:

\[ K_e = B^{\top}\,C\,B \cdot d \cdot A . \]

\(B\) übersetzt die sechs Eckverschiebungen in die drei Dehnungen — es sind die drei Summen von eben —, \(C\) die Dehnungen in Spannungen, und \(B^{\top}\) die Spannungen zurück in Knotenkräfte; \(d\) ist die Scheibendicke, \(A\) die Dreiecksfläche. Mehr steckt in der Formel nicht. Auch die Symmetrie, die Kapitel 9 angekündigt hat, sieht man ihr an: Lösung und Testfunktion stehen in \(\varepsilon^{\top} C\,\varepsilon\) gleichberechtigt, und was symmetrisch gebaut ist, kommt symmetrisch heraus. Die Farbkarte gleich unten zeigt es.

Welcher Knochen hier steckt. Ab jetzt wird mit echten Zahlen gerechnet, und die sind eine Entscheidung: Dieses und die folgenden Kapitel bis Kapitel 15 rechnen die Ferse als einen einzigen, homogen gedachten Knochen mit \(E = 17\) GPa und \(\nu = 0{,}30\) — den Werten der dichten Rinde, der Kortikalis (Kapitel 6). Der weiche Spongiosa-Kern im Inneren kommt darin nicht vor; das Modell ist also steifer als die echte Ferse, und seine Verschiebungen fallen zu klein aus. Wo im Knochen viel und wo wenig Spannung sitzt, hängt dagegen vor allem an Form und Kräften und bleibt lesbar. Welche Modelle dieses Buch sonst noch rechnet und wie sie sich unterscheiden, stellt Kapitel 13 nebeneinander.

Der folgende Block baut den Steckbrief für ein Dreieck des eben gezeigten Handnetzes (das Dreieck mit den Knoten 0, 9, 8) mit diesem Material und zeigt die 36 Einträge nicht als Zahlenwand, sondern als Abbildung 10.4: eine Farbkarte um null (rot positiv, blau negativ — das Vorzeichen ist die Richtung der Kraftantwort auf eine Einheitsverschiebung, kein Zug- oder Druckzustand), mit dem Zahlenwert klein in jeder Zelle — so werden die Symmetrie (Spiegelung an der Diagonale) und die verschwindenden Zeilensummen (die schmale Spalte rechts) auf einen Blick sichtbar. Zwei Eigenschaften rechnet der Block dabei ausdrücklich nach — die größte Unsymmetrie und die Zeilensumme über die waagerechten Freiheitsgrade; beide stehen darunter in Tabelle 10.2.

Den Steckbrief eines Handnetz-Dreiecks bauen und prüfen

# Ein Dreieck des Handnetzes, mit dem Material des Zahlenanhangs. Die Elementmatrix
# ist symmetrisch, und ihre Zeilensummen über die waagerechten bzw. senkrechten
# Freiheitsgrade sind null — eine reine Verschiebung des ganzen Dreiecks erzeugt
# keine inneren Kräfte (Starrkörperprobe).
dreieck = HANDNETZ_DREIECKE[1]                       # Knoten 0, 9, 8
ecken = [HANDNETZ_KNOTEN[dreieck[0]], HANDNETZ_KNOTEN[dreieck[1]],
         HANDNETZ_KNOTEN[dreieck[2]]]
Ke = elementmatrix(ecken, E_KORTIKAL, NU_KNOCHEN, TIEFE)
symmetrisch = max(abs(Ke[i][j] - Ke[j][i]) for i in range(6) for j in range(6))
zeilensumme_x = sum(Ke[0][j] for j in (0, 2, 4))     # waagerechte Beiträge, Zeile 0

# Die 36 Einträge als Farbkarte (in MN/m); Zeichenmechanik in programme/kap10.
zeichne_elementmatrix([[w / 1e6 for w in zeile] for zeile in Ke])

Abbildung 10.4: Der 6×6-Steckbrief eines Handnetz-Dreiecks (Knoten 0, 9, 8) als Farbkarte statt Zahlenwand. Zeile und Spalte sind die sechs Freiheitsgrade (u, v je Ecke); Eintrag (i, j) ist die Knotenkraft am Freiheitsgrad i bei Einheitsverschiebung von j. Divergente Skala um null (RdBu_r, symmetrisch in MN/m, an der Abbildung benannt): Rot positiv, Blau negativ. Die Karte ist an der Diagonale gespiegelt — die Matrix ist symmetrisch; die schmale Spalte rechts trägt je Zeile die Zeilensumme über die gleichgerichteten Freiheitsgrade und ist überall weiß (≈ 0) — eine Starrkörper-Translationsprobe (dass auch die starre Drehung keine inneren Kräfte erzeugt, prüft das Buch bei jedem Bau mit).

Tabelle 10.2: Die beiden eingebauten Kontrollen am Steckbrief des Dreiecks (0, 9, 8). Beide Werte sind Ausgaben der Rechnung, keine Behauptung.

Probe am Steckbrief	Wert (N/m)	Bedeutung
größte Unsymmetrie \(\|K_{ij} - K_{ji}\|\)	0,000000059605	Rundungsrest, also symmetrisch
Zeilensumme über die waagerechten Freiheitsgrade (Zeile 0)	0	null, also Starrkörperprobe bestanden

Deutung unter der Tabelle. Die Farbkarte zeigt es unmittelbar: Der Steckbrief ist symmetrisch — jede Zelle \((i,j)\) hat dieselbe Farbe und denselben Wert wie ihr Spiegelbild \((j,i)\); die größte gemessene Unsymmetrie liegt bei Rundungsresten. Das ist dieselbe Eigenschaft wie beim Feder-System aus Kapitel 4, und sie überträgt sich später auf die große Matrix \(K\). Die schmale Spalte rechts trägt die Zeilensummen über die gleichgerichteten Freiheitsgrade und ist durchweg weiß, also null: Schiebt man das ganze Dreieck starr zur Seite, ohne es zu verformen, entstehen keine inneren Kräfte. Das ist die eingebaute Kontrolle, dass der Steckbrief Physik und nicht Zufall enthält.

Ein Schnitt ist kein ganzer Knochen

In der Elementmatrix steckt eine Scheibendicke \(d\) — und darin eine bewusste Vereinfachung. Das ganze Buch rechnet nicht das dreidimensionale Fersenbein, sondern einen ebenen Sagittalschnitt (Seitenansicht) mit einer festen Ersatzdicke von 33 mm (Anhang A), und zwar unter der Annahme, die Kapitel 6 ebener Spannungszustand nennt: Die dünne Scheibe darf in die Tiefe frei atmen. Was quer zur Zeichenebene passiert — die Wölbung des Knochens, wie sich die Last über die echte Breite verteilt, die räumliche Anordnung der Trabekel — bildet das Modell nicht ab. Es beantwortet die Frage „wo im Schnitt wird es belastet?“ belastbar, aber es ist kein Ersatz für den ganzen Knochen. Wo diese Grenze wehtut, sagt das Buch es (etwa an der Verdrillung um die Längsachse, die ein ebener Schnitt grundsätzlich nicht sieht).

10.5 Der Zusammenbau: aus vielen kleinen Matrizen eine große

Jedes Dreieck liefert seinen \(6\times 6\)-Steckbrief. Die Assemblierung sortiert jeden seiner Einträge an die Stelle der beiden Freiheitsgrade ein, die er verbindet — und wo zwei Dreiecke sich einen Knoten teilen, addieren sich ihre Beiträge an derselben Stelle. Das lässt sich an den beiden hervorgehobenen Dreiecken aus Abbildung 10.3 Eintrag für Eintrag zeigen: das Dreieck \((0,9,8)\) und das Enthese-Dreieck \((9,10,8)\) — dessen Kante 9 → 10 die Enthese ist — teilen sich die Kante 8 → 9, also die beiden Knoten 8 und 9. Abbildung 10.5 stellt die beiden Dreiecke und ihre geteilte Kante noch einmal für sich heraus und zeigt daneben, welche Blöcke der Systemmatrix Beiträge aus beiden Dreiecken bekommen.

Code

import numpy as np
import matplotlib.pyplot as plt

kn = np.array(HANDNETZ_KNOTEN) * 1000.0
fig, (axl, axr) = plt.subplots(1, 2, figsize=(8.8, 4.2))

# Links: die zwei realen Dreiecke auf echten Koordinaten.
for dr, farbe in zip([(0, 9, 8), (9, 10, 8)], ("#e7c9a0", "#d9a066")):
    axl.fill(kn[list(dr), 0], kn[list(dr), 1], color=farbe, ec="#a9825f", lw=1.6)
axl.plot(kn[[8, 9], 0], kn[[8, 9], 1], color="#7a2b8a", lw=3.4,
         solid_capstyle="round")
for i in (0, 8, 9, 10):
    x, y = kn[i]
    axl.plot(x, y, "o", color="#5b6b7a", ms=14)
    axl.text(x, y, str(i), ha="center", va="center", fontsize=10, color="white")
orientierung(axl, mit_massstab=True)
axl.set_xlim(-6, 46); axl.set_ylim(-10, 34)
axl.set_aspect("equal"); axl.axis("off")

# Rechts: die 8×8-Blockstruktur über den vier Knoten; je Block ein Knotenpaar.
axr.set_title("Systemmatrix K (8×8)", fontsize=10)
knoten_reihe = [0, 8, 9, 10]

def n_beitraege(ni, nj):
    return sum(1 for dr in [(0, 9, 8), (9, 10, 8)] if ni in dr and nj in dr)

for bi, ni in enumerate(knoten_reihe):
    for bj, nj in enumerate(knoten_reihe):
        anz = n_beitraege(ni, nj)
        farbe = {0: "#ffffff", 1: "#f0e2d0", 2: "#e7d8ef"}[anz]
        rand = "#7a2b8a" if anz == 2 else "#cccccc"
        axr.add_patch(plt.Rectangle((bj, 3 - bi), 1, 1, fc=farbe, ec=rand,
                      lw=2.0 if anz == 2 else 0.6))
    axr.text(-0.3, 3 - bi + 0.5, str(ni), ha="right", va="center", fontsize=9,
             color="#555")
    axr.text(bi + 0.5, 4.12, str(ni), ha="center", va="bottom", fontsize=9,
             color="#555")
axr.text(2.0, -0.55, "violett: Beiträge beider Dreiecke", ha="center", fontsize=8,
         color="#7a2b8a")
axr.set_xlim(-0.9, 4.3); axr.set_ylim(-1.0, 4.6)
axr.set_aspect("equal"); axr.axis("off")
plt.tight_layout()
plt.show()

Abbildung 10.5: Assemblierung an den zwei konkreten Handnetz-Dreiecken (0, 9, 8) und (9, 10, 8). Links die beiden Dreiecke auf ihren echten Koordinaten; violett die geteilte Kante 8 → 9, deren beide Knoten 8 und 9 zu beiden Dreiecken gehören. Rechts die 8×8-Systemmatrix über den vier Knoten {0, 8, 9, 10}, in 2×2-Blöcke je Knotenpaar zerlegt: die Blöcke der geteilten Knoten 8 und 9 (violett) bekommen Beiträge aus BEIDEN Dreiecken addiert, die Blöcke von Knoten 0 (nur im ersten Dreieck) und Knoten 10 (nur im Enthese-Dreieck (9, 10, 8)) je nur aus einem (hellbraun); der Block, der die Knoten 0 und 10 koppeln würde, bleibt leer (weiß) — sie teilen sich kein Dreieck. Links Sagittalschnitt, x nach vorn, y nach oben; Maßstab 10 mm. Rechts schematisch.

Über die vier beteiligten Knoten \(\{0,8,9,10\}\) — acht Freiheitsgrade — bauen wir die kleine Systemmatrix jetzt ganz aus. Dazu genügt dieselbe assembliere-Funktion der oben; wir nummerieren die vier Knoten nur lokal um (0 → 0, 8 → 1, 9 → 2, 10 → 3), damit die \(8\times 8\)-Matrix genau diese vier Knoten trägt. Der folgende Block baut sie und stellt daneben die Frage, die den Zusammenbau ausmacht: Aus welchem Dreieck kommt eigentlich welcher Eintrag? Die Funktion beitraege_zu beantwortet sie, indem sie für eine Stelle der großen Matrix jedes Dreieck einzeln befragt — nur ein Dreieck, das beide Freiheitsgrade trägt, liefert überhaupt etwas.

Die zwei Dreiecke zusammenbauen — und nachsehen, woher jeder Eintrag kommt

# Die zwei benachbarten Dreiecke: (0, 9, 8) und das Enthese-Dreieck
# (9, 10, 8). Beide enthalten die Knoten 8 und 9 — die geteilte Kante. Über ihre
# vier Knoten {0, 8, 9, 10} bauen wir die kleine Systemmatrix; lokal umnummeriert
# 0 -> 0, 8 -> 1, 9 -> 2, 10 -> 3, damit die 8x8-Matrix genau diese vier trägt.
# Lokale Umnummerierung der vier Knoten (global -> lokal) und die zwei Dreiecke.
vier_knoten = [0, 8, 9, 10]
lokal = {}
for stelle, global_nr in enumerate(vier_knoten):
    lokal[global_nr] = stelle
lok_knoten = []
for global_nr in vier_knoten:
    lok_knoten.append(HANDNETZ_KNOTEN[global_nr])
lok_dreiecke = [(lokal[0], lokal[9], lokal[8]), (lokal[9], lokal[10], lokal[8])]
namen = ["u0", "v0", "u8", "v8", "u9", "v9", "u10", "v10"]

K8 = assembliere(lok_knoten, lok_dreiecke, E_KORTIKAL, NU_KNOCHEN, TIEFE)

def beitraege_zu(zeile, spalte):
    """Der Beitrag jedes der beiden Dreiecke zum 8×8-Eintrag [zeile][spalte]. Nur
    ein Dreieck, das BEIDE Freiheitsgrade trägt, liefert etwas — sonst bleibt der
    Eintrag null. So sieht man, welcher Eintrag aus welchem Dreieck kommt."""
    knoten_z = zeile // 2
    knoten_s = spalte // 2
    werte = []
    for dr in lok_dreiecke:
        if knoten_z in dr and knoten_s in dr:
            ecken = [lok_knoten[dr[0]], lok_knoten[dr[1]], lok_knoten[dr[2]]]
            Ke = elementmatrix(ecken, E_KORTIKAL, NU_KNOCHEN, TIEFE)
            pos_z = 2 * dr.index(knoten_z) + zeile % 2
            pos_s = 2 * dr.index(knoten_s) + spalte % 2
            werte.append(Ke[pos_z][pos_s])
    return werte

# Zwei Einträge des geteilten Blocks: die Diagonale des Knotens 8 und die
# Kopplung der beiden Endknoten 8 und 9 der geteilten Kante.
u8 = 2 * lokal[8]
u9 = 2 * lokal[9]
b_diag = beitraege_zu(u8, u8)
b_kante = beitraege_zu(u8, u9)

Alle 64 Einträge stehen in Tabelle 10.3, ihre Freiheitsgrade als Zeilen- und Spaltennamen (\(u,v\) je Knoten, Werte in MN/m); Tabelle 10.4 zerlegt danach zwei der geteilten Einträge in die Beiträge der beiden Dreiecke.

Tabelle 10.3: Die assemblierte 8×8-Systemmatrix der zwei Dreiecke (0, 9, 8) und (9, 10, 8) in MN/m. Die Zeilen und Spalten der geteilten Knoten 8 und 9 tragen Beiträge beider Dreiecke; die Kästchen, die die Knoten 0 und 10 koppeln würden, sind null.

	u0	v0	u8	v8	u9	v9	u10	v10
u0	307,9	73,5	−41,7	−92,5	−266,2	18,9	0	0
v0	73,5	146,2	−107,9	−119,2	34,3	−27	0	0
u8	−41,7	−107,9	282,6	−56,8	53,8	137,2	−294,7	27,5
v8	−92,5	−119,2	−56,8	428,1	137,2	−268,4	12,1	−40,5
u9	−266,2	34,3	53,8	137,2	628,5	25,8	−416,1	−197,3
v9	18,9	−27	137,2	−268,4	25,8	592,3	−181,9	−296,9
u10	0	0	−294,7	12,1	−416,1	−181,9	710,8	169,8
v10	0	0	27,5	−40,5	−197,3	−296,9	169,8	337,4

Tabelle 10.4: Zwei Einträge des geteilten Blocks, in ihre Herkunft zerlegt (MN/m). Die Spalte „Summe“ ist die Addition der beiden Dreiecksbeiträge, die Spalte „in der Matrix“ der Wert, den der Zusammenbau tatsächlich hingeschrieben hat — sie stimmen überein.

Eintrag	aus Dreieck (0, 9, 8)	aus Dreieck (9, 10, 8)	Summe	in der Matrix
\([u_8][u_8]\) — Diagonale des geteilten Knotens 8	113,6	169	282,6	282,6
\([u_8][u_9]\) — Kopplung der Kanten-Endknoten 8 und 9	−71,9	125,7	53,8	53,8

Deutung unter den Tabellen. Die \(8\times 8\)-Matrix trägt genau die vier Knoten der beiden Dreiecke. Der Block der geteilten Knoten 8 und 9 (die Zeilen und Spalten \(u_8, v_8, u_9, v_9\)) ist die Summe der Beiträge beider Dreiecke. Die erste Zeile von Tabelle 10.4 zeigt es an der Diagonale \([u_8][u_8]\) des geteilten Knotens 8; die zweite an der Kopplung \([u_8][u_9]\) der beiden Endknoten 8 und 9 der geteilten Kante — dieser Eintrag ist überhaupt nur besetzt, weil sich 8 und 9 eine Kante teilen, und auch er ist Zahl für Zahl die Summe beider Dreiecke. Die Blöcke von Knoten 0 (nur im Dreieck \((0,9,8)\)) und Knoten 10 (nur im Enthese-Dreieck \((9,10,8)\)) bekommen dagegen je nur einen Beitrag, und der Block, der die Knoten 0 und 10 koppeln würde, bleibt null — die beiden teilen sich kein Dreieck. Genau das ist der ganze Zusammenbau: einsortieren und an geteilten Knoten addieren.

Was hier an zwei Dreiecken geschieht, geschieht am ganzen Handnetz mit allen dreizehn: Jedes Dreieck legt seine 36 Einträge an die Zeilen und Spalten seiner drei Knoten, und wo Dreiecke sich Knoten teilen, summieren sich die Beiträge — so wächst die kleine \(8\times 8\) zur vollen \(26\times 26\)-Systemmatrix des Handnetzes. Das ist dieselbe Addition wie am Feder-System, wo zwei Federn am selben Knoten ihre Steifigkeiten addieren (Kapitel 4); und weil jeder Knoten nur mit seinen wenigen Netznachbarn koppelt, bleibt die große Matrix dünn besetzt. Wie das Netz sich Dreieck für Dreieck aufbaut und die Matrix sich dabei füllt, zeigt die folgende Animation.

Vorhersage-Punkt

Bevor du abspielst: Wenn beim Aufbau ein neues Dreieck an einen schon vorhandenen Knoten angesetzt wird — öffnet dieser geteilte Knoten neue Kästchen in der Matrix K, oder füllt das Dreieck nur die schon vorhandenen Zeilen und Spalten weiter? Lege dich fest.

Deutung. Das ist die Antwort auf die Vorhersagefrage: Ein geteilter Knoten tut beides. Auf seiner eigenen Zeile und Spalte (dem Diagonalkästchen) wird nur weiteraddiert — dort wächst der Wert. Zu den anderen Ecken des neuen Dreiecks aber öffnet er neue Kästchen, weil er jetzt mit diesen Knoten gekoppelt ist. Jedes Dreieck belegt so genau die Zellen seiner drei Knoten; die beiden inneren Knoten gehören zu den meisten Dreiecken (Knoten 8 zu neun, Knoten 11 zu sechs) und sammeln darum die meisten Einträge, die Randknoten koppeln nur mit ihren wenigen Nachbarn. So entsteht die dünn besetzte Systemmatrix, die der Löser aus Kapitel 5 knackt.

Wie dünn besetzt sie ist, lässt sich am Handnetz zählen. Das Besetzungsmuster einer Fersenmatrix ist in diesem Buch nicht neu — Abbildung 5.6 in Kapitel 5 hat es schon gezeigt, als Versprechen und noch ohne Erklärung, woher die Einträge kommen. Jetzt ist es bekannt: Jede besetzte Zelle stammt von mindestens einem Dreieck, das die beiden Knoten dieser Zeile und Spalte gemeinsam hat. Am Handnetz mit seinen dreizehn Knoten kann man das Kästchen für Kästchen nachsehen, und der folgende Block tut genau das — er baut mit derselben assembliere-Funktion die volle \(26\times 26\)-Matrix des ganzen Handnetzes und zeichnet nur noch, welche Zelle überhaupt einen Eintrag bekommt. Zum Vergleich zählt er dieselbe Größe am groben Buchnetz mit; Tabelle 10.5 hält beides nebeneinander.

Die volle Handnetz-Matrix bauen und ihr Besetzungsmuster zeigen

# Alle dreizehn Dreiecke, dieselbe assembliere-Funktion wie oben: die volle
# Systemmatrix des Handnetzes. Gezeigt wird nur, WELCHE Zelle einen Eintrag
# bekommt — nicht, welchen Wert.
K_hand = assembliere(HANDNETZ_KNOTEN, HANDNETZ_DREIECKE,
                     E_KORTIKAL, NU_KNOCHEN, TIEFE)
besetzt_hand = np.array([[abs(wert) > 0.0 for wert in zeile] for zeile in K_hand])
fg_hand = len(K_hand)
fuell_hand = 100.0 * besetzt_hand.sum() / (fg_hand * fg_hand)

# Dieselbe Zählung am groben Buchnetz — nur die Zahlen, ohne Bild.
netz = BUCHNETZE["grob"]
_, fuell_grob = besetzung_gross(netz["dreiecke"], netz["n_knoten"])
fg_grob = 2 * netz["n_knoten"]

fig, ax = plt.subplots(figsize=(5.8, 5.8))
ax.imshow(~besetzt_hand, cmap="gray", vmin=0, vmax=1, interpolation="nearest")
for k in range(1, len(HANDNETZ_KNOTEN)):
    ax.axhline(2 * k - 0.5, color="#c9a37a", lw=0.7)
    ax.axvline(2 * k - 0.5, color="#c9a37a", lw=0.7)
mitten = [2 * k + 0.5 for k in range(len(HANDNETZ_KNOTEN))]
namen_kn = [str(k) for k in range(len(HANDNETZ_KNOTEN))]
ax.set_xticks(mitten); ax.set_xticklabels(namen_kn, fontsize=8)
ax.set_yticks(mitten); ax.set_yticklabels(namen_kn, fontsize=8)
ax.set_xlabel("Knoten (je zwei Freiheitsgrade)", fontsize=9)
ax.set_ylabel("Knoten (je zwei Freiheitsgrade)", fontsize=9)
ax.set_title(f"Handnetz: {fg_hand}×{fg_hand}-Matrix", fontsize=10, pad=10)
plt.tight_layout()
plt.show()

Abbildung 10.6: Besetzungsmuster der vollen Systemmatrix K des Handnetzes (13 Knoten, 26 Freiheitsgrade): dunkel jede Zelle, die einen Eintrag bekommt, weiß der Rest. Die hellen Linien trennen die 2×2-Blöcke der Knotenpaare, die Zahlen an den Rändern sind die Knotennummern aus Abbildung 10.3. Ein Block ist genau dann besetzt, wenn seine beiden Knoten sich mindestens ein Dreieck teilen: Die Zeile des inneren Knotens 8, der zu neun der dreizehn Dreiecke gehört, ist fast voll; die Zeile von Knoten 12, der nur an 6, 7 und 8 hängt, hat vier besetzte Blöcke — sich selbst und diese drei. Kein Feldbild, sondern die Struktur der Matrix.

Tabelle 10.5: Zwei Netze, dieselbe Bauart, sehr verschiedene Füllgrade. Alles, was nicht besetzt ist, ist exakt null — der Löser muss es nie speichern und nie anfassen.

Systemmatrix	Freiheitsgrade	Plätze	besetzt
Handnetz (13 Knoten)	26	676	37,3 %
grobes Buchnetz (330 Knoten)	660	0,44 Millionen	1,98 %

Deutung unter der Tabelle. Am Handnetz ist noch gut ein Drittel der Matrix besetzt — bei dreizehn Knoten ist eben fast jeder mit fast jedem über irgendein Dreieck verbunden. Am groben Buchnetz mit seinen 330 Knoten sind es nur noch rund zwei Prozent. Das ist der Punkt: Ein Knoten koppelt immer nur mit seinen wenigen direkten Netznachbarn, und diese Zahl bleibt beim Verfeinern ungefähr gleich, während die Zahl der Plätze mit dem Quadrat der Knotenzahl wächst. Je größer das Netz, desto leerer die Matrix — genau die Eigenschaft, die Kapitel 5 angekündigt hat und die große FEM-Systeme überhaupt erst handhabbar macht.

Und die rechte Seite? \(K\) ist nur die halbe Miete; das Gleichungssystem \(K\,u = f\) braucht auch ein \(f\). Es entsteht auf demselben Weg. In der schwachen Form von Kapitel 9 stand rechts \(\int A\,b\,v\,\mathrm{d}x\) — die Last, gewichtet mit der Testfunktion. In der Ebene wird daraus ein Integral über die Dreiecke (für eine Last im Inneren) und über Randkanten (für eine Last, die von außen angreift), und beides wird genauso einsortiert wie die Steifigkeiten: Jeder Beitrag landet an den Freiheitsgraden der Knoten, die er berührt, und wo mehrere Beiträge denselben Knoten treffen, addieren sie sich. Das Fersenmodell dieses Buches kennt keine Last im Inneren — das Eigengewicht des Knochens ist gegen die Schrittkräfte verschwindend —, sein \(f\) kommt also allein von den Rändern. Kapitel 11 setzt zum ersten Mal eine solche Randlast ein, Kapitel 12 behandelt sie systematisch.

10.6 Dreiecksform am Regler: was spitze Dreiecke anrichten

Nicht nur die Zahl der Dreiecke zählt, auch ihre Form. Der Regler nimmt einen winzigen Patch aus vier Dreiecken um einen frei beweglichen Innenknoten und schiebt diesen zum Rand, bis die Dreiecke zu spitzen Splittern werden — und rechnet dabei live die Formgüte, den größten Elementmatrix-Eintrag, die Konditionszahl der reduzierten Systemmatrix und die Senkung des Knotens.

Vorhersage-Punkt

Bevor du schiebst: Sehr spitze, schlecht geformte Dreiecke — bleibt die Rechnung stabil, oder wird das Ergebnis unzuverlässig? Lege dich fest.

Deutung. Je spitzer die Dreiecke, desto kleiner die Formgüte (von rund \(0{,}87\) für die anfangs rechtwinklig-gleichschenkligen Dreiecke gegen null) und desto größer der größte Eintrag der Elementmatrix — er wächst umgekehrt zur Fläche, also ohne obere Schranke. Die Antwort auf die Vorhersagefrage ist zweigeteilt und genau das Lehrreiche: Dieses winzige System bleibt stabil — seine Konditionszahl steigt nur von rund eins auf rund zwei, die Rechnung löst sich weiter gutmütig. Die Gefahr zeigt sich also nicht an einem einzelnen spitzen Dreieck, sondern erst im großen Zusammenbau: Dort mischen sich die über alle Grenzen wachsenden Einträge der Splitter mit den winzigen der gut geformten Dreiecke in derselben Matrix, und dieses Mischen von Riesig und Winzig verschlechtert die Kondition der Gesamtmatrix — bis die Lösung unzuverlässig wird. Deshalb erzeugen die Netzgeneratoren bewusst gut geformte Dreiecke, und deshalb ist die Form ein Qualitätsmerkmal des Netzes, kein Schönheitsfehler.

Editierbare Zelle. Der Regler zeigt, dass der größte Eintrag wächst. Offen bleibt, wie teuer eine bestimmte Formgüte wird — und das lässt sich ausrechnen. Der folgende Block formt dasselbe Dreieck über eine Reihe von Spitzigkeiten, sucht die Stelle, an der die Formgüte unter eine selbst gewählte Schranke fällt, und sagt, um welches Vielfache der größte Elementmatrix-Eintrag dort gegenüber dem gut geformten Dreieck gewachsen ist. Verändere guete_schwelle und führe die Zelle aus: Bei welcher Schranke ist der Aufpreis noch erträglich?

Deutung. Die beiden Kurven laufen gegeneinander: Die blaue Formgüte fällt stetig gegen null, während der rote Faktor erst gemächlich, dann immer steiler steigt. Die gestrichelte Senkrechte markiert die gewählte Schranke, die Ausgabe darüber nennt Spitzigkeit und Faktor als Zahlen. Mit der voreingestellten Schranke \(0{,}3\) ist der Aufpreis noch bescheiden — knapp das Vierfache. Erst dicht an der Entartung wird es teuer: bei Güte \(0{,}1\) schon rund das Elffache, und weiter ohne obere Schranke, weil der größte Eintrag umgekehrt zur Fläche wächst. Genau deshalb geben Netzgeneratoren eine Mindestgüte vor, statt jedes spitze Dreieck zu verbieten: Ein bisschen Spitze kostet wenig, ein Splitter kostet alles.

Übungen

Ü 10.1 (Verstehen). Zeichne zwei benachbarte Dreiecke, die sich einen Knoten teilen, und markiere für diesen Knoten die Zelle in der großen Matrix \(K\). Erkläre in zwei, drei Sätzen, warum genau diese Zelle die Summe der Beiträge beider Dreiecke trägt und warum \(K\) deshalb dünn besetzt ist.

Musterlösung zu Ü 10.1

Jeder Eintrag der Elementmatrix gehört zu einem Paar von Freiheitsgraden. Der geteilte Knoten hat in beiden Dreiecken denselben globalen Freiheitsgrad, also dieselbe Zeile und Spalte in \(K\). Beim Zusammenbau werden die beiden Beiträge an dieser Stelle addiert — wie zwei Federn, die am selben Knoten ziehen (Kapitel 4). \(K\) ist dünn besetzt, weil ein Knoten nur mit den wenigen Knoten koppelt, mit denen er sich ein Dreieck teilt; zu allen anderen bleibt die Zelle null.

Ü 10.2 (Verändern). Öffne den Dreiecksform-Regler und protokolliere Formgüte, größten Elementmatrix-Eintrag und die Kondition der reduzierten Matrix bei Spitzigkeit 0, 0,5 und 0,95. Sage vorher, ob das kleine System dabei instabil wird. Was beobachtest du — und warum warnt das Kapitel trotzdem vor spitzen Dreiecken?

Musterlösung zu Ü 10.2

Je spitzer, desto kleiner die Güte und desto größer der größte Elementmatrix-Eintrag — er wächst umgekehrt zur Fläche, ohne obere Schranke. Die Kondition dieses winzigen 2×2-Systems steigt aber nur mäßig, von rund 1 auf rund 2; die Rechnung bleibt hier gutmütig. Die Gefahr zeigt sich also nicht am einzelnen spitzen Dreieck, sondern erst im großen Zusammenbau, wo sich die über alle Grenzen wachsenden Splitter-Einträge mit den winzigen der gut geformten Dreiecke in derselben Matrix mischen und die Kondition der Gesamtmatrix verschlechtern. Deshalb erzeugen die Netzgeneratoren bewusst gut geformte Dreiecke.

Ü 10.3 (Übertragen). Begründe, warum das Handnetz zum Verstehen taugt, aber nicht für die Zahlen des Buches. Was muss ein Netz können, damit man seinen Zahlen traut?

Musterlösung zu Ü 10.3

Das Handnetz ist so klein, dass man jeden Schritt von Hand nachrechnen kann — Elementmatrix, Zusammenbau und, in Kapitel 11, auch das Lösen unter Last. Genau dafür ist es da. Aber dreizehn Knoten lösen die Spannung an einer Lasteinleitung viel zu grob auf: Dort ändert sich die Belastung über wenige Millimeter stark, und ein Dreieck kennt nur eine einzige Dehnung. Die Zahlen des Buches stammen deshalb von den feinen Netzen hinter den Reglern und aus einer getrennt gerechneten Referenzrechnung. Dass das Handnetz und der schnelle Rechner hinter den Reglern auf demselben kleinen Netz dieselben Zahlen liefern, prüft das Buch bei jedem Bau — das sichert die Methode, nicht die Feinheit. Wie fein ein Netz sein muss, ist die Frage von Kapitel 13.

Roter Faden

Wo kam das schon vor, wo kommt es wieder? Zurück: Die schwache Form aus Kapitel 9 wird hier konkret — die stückweise linearen Ansätze sitzen jetzt auf Dreiecken, und weil sie dort konstante Dehnungen liefern, schrumpft ihr Integral zu „Wert mal Fläche”: der Steckbrief \(K_e = B^{\top} C B \cdot d \cdot A\). Die Materialtabelle \(C\) darin ist die 3×3-Tabelle aus Kapitel 6, mit \(E\) und \(\nu\) von dort; die Symmetrie des Steckbriefs ist die Galerkin-Symmetrie aus Kapitel 9. Die Elementmatrizen bauen dieselbe Steifigkeitsmatrix \(K\) wie das Federsystem aus Kapitel 4, gelöst mit dem Gauß-Verfahren aus Kapitel 5; das Besetzungsmuster, das Kapitel 5 vorweggenommen hat, ist jetzt Kästchen für Kästchen erklärt. In der Studiums-Tabelle im Vorwort steht dafür die Zeile „Netz aus Dreiecken → Diskretisierung, Ansatzfunktionen” — mit diesem Kapitel ist sie eingelöst; fortgeschrieben wird die Tabelle am Ende von Teil IV (Kapitel 13). Vor: Kapitel 11 rechnet dieses Handnetz erstmals vollständig durch — aus derselben Assemblierung, dem Lastvektor eines Schrittmoments aus dem Zahlenanhang und dem Talus-Lager entsteht das erste Feldbild aus reiner Handrechnung. Kapitel 12 führt die vier Randbedingungen (Boden, Sehne, Faszie, Lager) danach systematisch am feinen Netz ein; Kapitel 13 prüft, wie fein das Netz sein muss. Das Netz ist die Bühne für alles Weitere.

Was dieses Kapitel NICHT tut

Es baut keine Elemente höherer Ordnung (nur P1-Dreiecke) und keine gekrümmten Ränder; es vernetzt nicht im Browser (die Netze erzeugt gmsh offline) und verfeinert nicht adaptiv. Es löst das Handnetz noch nicht unter einer echten Last — die vollständige Rechnung mit Lastvektor und Talus-Lager steht in Kapitel 11. Es hängt keine echten Randbedingungen ans Netz (Boden, Sehne, Faszie, Talus-Lager kommen in Kapitel 12) und klärt nicht, wie fein das Netz sein muss, noch rechnet es eine Fehleranalyse — das ist Kapitel 13. Der Dreiecksform-Patch und das einzelne Dreieck rechnen ohne Einheit und ohne Bezug zur Ferse; echte Geometrie und echtes Material tragen nur das Handnetz und die drei Buchnetze.