10 In die Fläche: Dreiecke

Ein Mosaik pflastert jede Wand — krumme Bögen, Ecken, Löcher —, indem es kleine Steinchen dicht an dicht legt. Ein Stadiondach, eine geodätische Kuppel, das Drahtgittermodell eines Kinofilmgesichts: Schaut man genau hin, sind sie aus Dreiecken gebaut. Das ist kein Zufall, sondern der Grund, warum Dreiecke das Lieblingswerkzeug der Finite-Elemente-Methode sind: Mit genug Dreiecken lässt sich jede Form auslegen.

Bis hierher lebte unser Träger auf einer Linie — fünf Knoten, vier Elemente, eine Zahl pro Ort. In diesem Kapitel verlässt er die Linie und wird zur Fläche. Und das Beste kommt gleich vorweg: Es ist dasselbe Spiel. Alles, was du seit Kapitel 3 gelernt hast — Steckbrief, Assemblieren, Randwerte, Lösen —, funktioniert unverändert weiter. Es ändert sich nur eine einzige Frage: Wie interpoliert man auf einem Dreieck?

Code

import matplotlib.pyplot as plt
import numpy as np

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10.5, 3.8))

# Links: Viertelkreis mit Dreiecken ausgelegt
ringe = 4
punkte = []
for r in range(ringe + 1):
    radius = r / ringe
    for w in range(r + 1):
        winkel = (np.pi / 2) * (w / r) if r > 0 else 0.0
        punkte.append((radius * np.cos(winkel), radius * np.sin(winkel)))
xs = [p[0] for p in punkte]
ys = [p[1] for p in punkte]
ax1.triplot(xs, ys, color="tab:blue", lw=0.8)
kreis = np.linspace(0, np.pi / 2, 60)
ax1.plot(np.cos(kreis), np.sin(kreis), color="tab:orange", lw=2.0)
ax1.set_title("Jede Form lässt sich triangulieren", fontsize=10)
ax1.set_aspect("equal")
ax1.axis("off")

# Rechts: der Träger mit 8 Dreiecken
knoten = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0),
          (0.0, 0.05), (0.5, 0.05), (1.0, 0.05),
          (0.0, 0.10), (0.5, 0.10), (1.0, 0.10)]
dreiecke = [(0, 1, 4), (0, 4, 3), (1, 2, 5), (1, 5, 4),
            (3, 4, 7), (3, 7, 6), (4, 5, 8), (4, 8, 7)]
bx = [p[0] for p in knoten]
by = [p[1] for p in knoten]
ax2.triplot(bx, by, dreiecke, color="tab:blue", lw=1.0)
ax2.plot(bx, by, "o", color="tab:orange", ms=8, markeredgecolor="black")
ax2.set_title("Unser Träger, mit 8 Dreiecken gepflastert", fontsize=10)
ax2.set_xlabel("Länge x (m)")
ax2.set_aspect("equal")
plt.tight_layout()
plt.show()

Abbildung 10.1: Warum Dreiecke. Links eine krumme Form (ein Viertelkreis), mit Dreiecken ausgelegt — die Ränder werden treppchenweise erfasst, mit mehr Dreiecken beliebig genau. Rechts derselbe Gedanke für unseren Träger: das grobe 8-Dreiecke-Netz dieses Kapitels. Drei Punkte legen immer eine Ebene fest, darum ist das Dreieck der einfachste Flächen-Baustein, den es gibt.

Was man hier sieht: links eine gekrümmte Form, aus lauter Dreiecken gelegt — je feiner die Dreiecke, desto glatter der Rand. Rechts der Träger dieses Kapitels: neun Knoten, acht Dreiecke, klein genug, dass wir sein Gleichungssystem gleich noch von Hand anschauen. Beide Male gilt: Drei Punkte spannen genau eine Ebene auf, und über einer Ebene lässt sich kinderleicht interpolieren. Das ist der ganze Trick.

Lernziele

Nach diesem Kapitel kannst du …

… erklären, warum in der Fläche Dreiecke die natürlichen Bausteine sind — jede Form lässt sich triangulieren, und drei Ecken legen eine Ebene fest,
… einen Punkt im Dreieck durch seine drei Flächenanteile beschreiben und damit zwischen drei Eckwerten interpolieren,
… die Ansatzfunktion eines Knotens als Zelt über seinen Dreiecken beschreiben — die flächenhafte Fortsetzung des Hütchens aus Kapitel 7,
… ein Netz als Knotenliste plus Dreiecksliste lesen und die 3×3-Elementmatrix eines Dreiecks aus den Steigungen seiner Zeltflächen aufstellen,
… am 8-Dreiecke-Träger die Assemblierung nachvollziehen, das System mit dem Löser aus Kapitel 4 knacken und die Bandstruktur deuten.

Naheliegende Vermutung

Vermutung: „Zwei Dimensionen sind grundsätzlich komplizierter als eine. Für die Fläche braucht man ganz neue Mathematik — die Werkzeuge aus dem 1D-Teil kann ich vergessen.”

Warum sie naheliegt: Eine Fläche hat mehr Richtungen, mehr Nachbarn, mehr zu bedenken als eine Linie. Es fühlt sich nach einem großen Sprung an, und in vielen Lehrbüchern beginnt hier ein neues, dickeres Kapitel voller Indizes.

Was stattdessen stimmt: Die Rezeptur ist Wort für Wort dieselbe wie in 1D — Ansatz, Elementmatrix, Assemblierung, Randwerte, Lösen. Es kommt genau ein neuer Baustein hinzu: Interpolation mit Flächenanteilen statt Streckenanteilen. Eine Idee, nicht zehn. Wer Kapitel 7 verstanden hat, hat den Berg schon bestiegen; dieses Kapitel ist der Spaziergang auf dem Kamm.

10.1 Der Träger als Fläche: die Geometrie festlegen

Bevor wir rechnen, klären wir, was die Fläche eigentlich zeigt. Unser Träger ist ein Quader: 1 m lang, 0,1 m hoch, 0,1 m tief. In 1D haben wir ihn auf seine Länge zusammengedrückt und Höhe und Tiefe in der Querschnittsfläche \(A = 0{,}01\ \mathrm{m^2}\) versteckt. Jetzt lösen wir eine der beiden versteckten Richtungen auf: Wir schauen von der Seite auf den Träger.

Code

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(9.0, 3.4))

# Der Träger von der Seite (x-y-Ebene)
ax.add_patch(plt.Rectangle((0, 0), 1.0, 0.1, facecolor="#fde8d0",
                           edgecolor="black", lw=1.5))
# Tiefen-Pfeil (in die Seite hinein, schräg gezeichnet)
ax.annotate("", xy=(1.12, 0.14), xytext=(1.0, 0.02),
            arrowprops=dict(arrowstyle="-|>", color="0.4", lw=2.0))
ax.text(1.14, 0.15, "Tiefe d = 0,1 m\n(in die Seite hinein)", color="0.4",
        fontsize=9, va="bottom")

# Massketten
ax.annotate("", xy=(1.0, -0.03), xytext=(0.0, -0.03),
            arrowprops=dict(arrowstyle="<|-|>", color="tab:blue", lw=1.3))
ax.text(0.5, -0.055, "Länge x: 0 … 1 m", color="tab:blue", ha="center",
        fontsize=10)
ax.annotate("", xy=(-0.03, 0.1), xytext=(-0.03, 0.0),
            arrowprops=dict(arrowstyle="<|-|>", color="tab:green", lw=1.3))
ax.text(-0.05, 0.05, "Höhe y:\n0 … 0,1 m", color="tab:green", ha="right",
        va="center", fontsize=10)

# kalt links, Brandherd rechts
ax.add_patch(plt.Rectangle((-0.11, 0.0), 0.03, 0.1, hatch="///",
                           facecolor="lightsteelblue", edgecolor="black"))
ax.text(-0.095, 0.12, "20 °C", color="tab:blue", ha="center", fontsize=9)
ax.fill([1.0, 1.06, 1.0], [0.03, 0.05, 0.07], color="orangered")
ax.text(1.03, 0.115, "Brandherd\n300 °C", color="orangered", ha="center",
        fontsize=9)

ax.set_xlim(-0.25, 1.35)
ax.set_ylim(-0.12, 0.32)
ax.set_aspect("equal")
ax.axis("off")
plt.tight_layout()
plt.show()

Abbildung 10.2: Die Geometrie-Konvention des ganzen restlichen Buches. Wir sehen den Träger von der Seite: x ist die Länge (0 bis 1 m), y die Höhe (0 bis 0,1 m). Die dritte Richtung, die Tiefe d = 0,1 m, steht senkrecht auf dem Papier (grauer Pfeil). Höhe mal Tiefe ergibt weiterhin die 1D-Querschnittsfläche A = 0,1 m · 0,1 m = 0,01 m².

Was man hier sieht: den Träger von der Seite. Die Länge \(x\) läuft von links (kühles Auflager, 20 °C) nach rechts (Brandherd, 300 °C), die Höhe \(y\) von unten nach oben. Die Tiefe \(d = 0{,}1\ \mathrm{m}\) zeigt in die Papierebene hinein — sie bleibt in diesem Buch immer gleich. Damit stimmt alles mit dem 1D-Teil zusammen: Höhe mal Tiefe ist dieselbe Querschnittsfläche wie in 1D, \(0{,}1 \cdot 0{,}1 = 0{,}01\ \mathrm{m^2} = A\).

Merke: die Geometrie-Konvention

Das 2D-Modell ist die Seitenansicht des Trägers:

\(x\) = Länge, von 0 bis 1 m,
\(y\) = Höhe, von 0 bis 0,1 m,
Tiefe \(d\) = 0,1 m (senkrecht zur Zeichenebene, immer konstant).

Wo in 1D die ganze Querschnittsfläche \(A = 0{,}01\ \mathrm{m^2}\) als ein Faktor auftrat, tritt in 2D nur noch die Tiefe \(d = 0{,}1\ \mathrm{m}\) auf — die Höhe ist jetzt eine echte Richtung des Netzes und keine versteckte Zahl mehr. Diese Konvention gilt bis zum Finale (Kapitel 14, Kapitel 15).

Jetzt, wo die Höhe eine echte Richtung des Netzes wird, dürfen wir eine Frage einlösen, die im 1D-Teil offengeblieben war: Warum durften wir den Träger überhaupt als bloße Linie behandeln? In Kapitel 9 stand die Antwort schon als Vorgriff bereit — hier ist sie, und mit dem 2D-Modell können wir sie gleich nachprüfen.

Warum durften wir 1D rechnen?

Der ganze dritte Teil hat den Träger als Linie gerechnet — eine Temperatur je Ort entlang der Länge, die Höhe zu einer bloßen Zahl in der Querschnittsfläche eingedampft. Das war kein Trick, sondern durfte sein, und der Grund steckt in den Zeitkonstanten aus Kapitel 9.

Ein Temperaturunterschied gleicht sich über eine Strecke \(\ell\) in etwa der Zeit \(\tau = \ell^2\rho c/\lambda\) aus — die Strecke geht im Quadrat ein. Über die volle Länge \(L = 1\ \mathrm{m}\) des Trägers ist das

\[ \tau_{\text{längs}} = \frac{L^2\,\rho\,c}{\lambda} = \frac{1^2 \cdot 7850 \cdot 460}{50} \approx 72\,220\ \mathrm{s} \approx 20\ \mathrm{h} , \]

über die kleine Höhe \(H = 0{,}1\ \mathrm{m}\) dagegen

\[ \tau_{\text{Höhe}} = \frac{H^2\,\rho\,c}{\lambda} = \frac{0{,}1^2 \cdot 7850 \cdot 460}{50} \approx 722\ \mathrm{s} \approx 12\ \mathrm{min} . \]

Das Verhältnis ist genau

\[ \frac{\tau_{\text{längs}}}{\tau_{\text{Höhe}}} = \left(\frac{L}{H}\right)^2 = \left(\frac{1}{0{,}1}\right)^2 = 100 : \]

Quer über die Höhe gleicht sich die Temperatur rund hundertmal schneller aus als längs. Ehe sich über die Länge etwas rührt, ist der Träger über seine Höhe längst durchgewärmt — auf jedem senkrechten Schnitt herrscht praktisch eine einzige Temperatur. Genau das rechtfertigt die 1D-Sicht: eine Zahl je Längsstelle genügt.

Und jetzt, mit dem Netz in der Fläche, müssen wir es nicht mehr glauben — wir können es nachprüfen. Gleich wird die Rechnung mit gleichem Rand links und rechts von selbst ein Feld liefern, das über die Höhe keinen Unterschied macht. Erst wenn die Brandgase die Höhe ungleich belasten — von unten heiß, von oben kühl —, wird der Höhenunterschied wichtig; das ist die Szene, die den Träger im Finale krümmt (Kapitel 14).

10.2 Vom Streckenanteil zum Flächenanteil

Erinnern wir uns, wie Interpolation in 1D funktioniert hat. Zwischen zwei Knoten teilt ein Punkt die Strecke in zwei Anteile: Liegt er bei einem Viertel des Weges, gehört er zu drei Vierteln dem linken und zu einem Viertel dem rechten Knoten. Genau das haben die Hütchen aus Kapitel 7 aufgeschrieben — der Anteil eines Knotens ist das Gewicht seines Wertes.

In der Fläche ist ein Punkt nicht mehr durch eine Zahl (den Streckenanteil) festgelegt, sondern von drei Ecken umgeben. Der Trick: Verbinde den Punkt \(P\) mit den drei Ecken. Das zerlegt das Dreieck in drei Teildreiecke, und die drei Flächenanteile übernehmen exakt die Rolle der Streckenanteile.

Code

import matplotlib.pyplot as plt

V1 = (0.0, 0.0)
V2 = (4.0, 0.0)
V3 = (0.0, 4.0)
P = (1.0, 1.0)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10.5, 4.2),
                               gridspec_kw={"width_ratios": [1.5, 1]})

# Teildreiecke einfärben
ax1.fill([P[0], V2[0], V3[0]], [P[1], V2[1], V3[1]], color="tab:blue",
         alpha=0.45)
ax1.fill([P[0], V3[0], V1[0]], [P[1], V3[1], V1[1]], color="tab:orange",
         alpha=0.45)
ax1.fill([P[0], V1[0], V2[0]], [P[1], V1[1], V2[1]], color="tab:green",
         alpha=0.45)
# Verbindungslinien und Rand
for V in (V1, V2, V3):
    ax1.plot([P[0], V[0]], [P[1], V[1]], color="0.4", lw=1.0)
ax1.plot([V1[0], V2[0], V3[0], V1[0]], [V1[1], V2[1], V3[1], V1[1]],
         color="black", lw=1.8)
# Beschriftung
ax1.text(1.6, 1.6, "A₁ = 4", color="tab:blue", fontsize=11)
ax1.text(0.15, 1.3, "A₂\n= 2", color="tab:orange", fontsize=11)
ax1.text(1.0, 0.18, "A₃ = 2", color="tab:green", fontsize=11)
ax1.plot(*P, "o", color="black", ms=9, zorder=5)
ax1.text(P[0] + 0.12, P[1] + 0.12, "P = (1, 1)", fontsize=10)
ax1.text(V1[0] - 0.15, V1[1] - 0.35, "Ecke 1\n(0, 0)", ha="center", fontsize=9)
ax1.text(V2[0], V2[1] - 0.35, "Ecke 2\n(4, 0)", ha="center", fontsize=9)
ax1.text(V3[0] - 0.15, V3[1] + 0.1, "Ecke 3 (0, 4)", ha="left", fontsize=9)
ax1.set_xlim(-0.7, 4.5)
ax1.set_ylim(-0.7, 4.5)
ax1.set_aspect("equal")
ax1.set_title("Punkt P teilt das Dreieck in drei Teilflächen", fontsize=10)
ax1.axis("off")

# Balkendiagramm der Anteile
anteile = [0.5, 0.25, 0.25]
farben = ["tab:blue", "tab:orange", "tab:green"]
namen = ["Anteil 1", "Anteil 2", "Anteil 3"]
ax2.bar(namen, anteile, color=farben)
for i, a in enumerate(anteile):
    ax2.text(i, a + 0.02, "%.0f %%" % (a * 100), ha="center", fontsize=10)
ax2.set_ylim(0, 0.65)
ax2.set_ylabel("Flächenanteil")
ax2.set_title("… und das sind die Gewichte", fontsize=10)
plt.tight_layout()
plt.show()

Abbildung 10.3: Der einzige neue Baustein des Kapitels. Verbindet man den Punkt P mit den drei Ecken, zerfällt das Dreieck in drei farbige Teildreiecke. Der Anteil einer Ecke ist die Fläche des Teildreiecks GEGENÜBER von ihr, geteilt durch die Gesamtfläche. Rechts die drei Anteile als Balken — sie summieren sich immer zu 100 %.

Was man hier sieht: Der Punkt \(P\) zerteilt das Dreieck in drei Teildreiecke, jedes in einer Farbe. Die Überkreuz-Logik ist wichtig: Der Anteil, der zu Ecke 1 gehört, ist die Fläche des Teildreiecks gegenüber von Ecke 1 (das blaue, das die Ecke 1 gar nicht berührt). Warum gegenüber? Wenn \(P\) direkt auf Ecke 1 säße, wäre das gegenüberliegende Teildreieck das ganze Dreieck — Anteil 1 wäre 100 %. Wandert \(P\) weg von Ecke 1, schrumpft dieses Teildreieck, und der Anteil sinkt.

10.2.1 Rechnen mit echten Zahlen

Die Flächen berechnen wir mit einem Schulrezept, das ganz ohne Vektorprodukt auskommt — der halben Kreuzdifferenz. Für ein Dreieck mit den Ecken \((x_1, y_1)\), \((x_2, y_2)\), \((x_3, y_3)\):

\[ \text{Fläche} = \tfrac{1}{2}\,\bigl|\, (x_2 - x_1)(y_3 - y_1) - (x_3 - x_1)(y_2 - y_1)\,\bigr| . \]

In Worten: Spanne von der ersten Ecke aus die beiden Kanten auf, bilde „erste Kante breit mal zweite Kante hoch, minus über Kreuz”, und nimm die Hälfte des Betrags. Für unser Dreieck mit \(V_1 = (0,0)\), \(V_2 = (4,0)\), \(V_3 = (0,4)\):

\[ A = \tfrac{1}{2}\,\bigl|(4-0)(4-0) - (0-0)(0-0)\bigr| = \tfrac{1}{2}\cdot 16 = 8 . \]

Jetzt die drei Teildreiecke mit \(P = (1,1)\), jedes mit derselben Formel:

\[ \begin{aligned} A_1 &= \text{Fläche}(P, V_2, V_3) = \tfrac{1}{2}\bigl|(4-1)(4-1) - (0-1)(0-1)\bigr| = \tfrac{1}{2}\,|9 - 1| = 4,\\[2pt] A_2 &= \text{Fläche}(P, V_3, V_1) = \tfrac{1}{2}\bigl|(0-1)(0-1) - (0-1)(4-1)\bigr| = \tfrac{1}{2}\,|1 + 3| = 2,\\[2pt] A_3 &= \text{Fläche}(P, V_1, V_2) = \tfrac{1}{2}\bigl|(0-1)(0-1) - (4-1)(0-1)\bigr| = \tfrac{1}{2}\,|1 + 3| = 2 . \end{aligned} \]

Die drei Flächenanteile sind die Teilflächen, geteilt durch die Gesamtfläche:

\[ \lambda_1 = \frac{A_1}{A} = \frac{4}{8} = 0{,}5, \qquad \lambda_2 = \frac{A_2}{A} = \frac{2}{8} = 0{,}25, \qquad \lambda_3 = \frac{A_3}{A} = \frac{2}{8} = 0{,}25 . \]

Die Probe ist immer dieselbe: \(0{,}5 + 0{,}25 + 0{,}25 = 1\) — die Anteile summieren sich zu 100 %. (Ich schreibe die Anteile hier \(\lambda_1, \lambda_2, \lambda_3\); das ist nur ein Name, kein neuer Stoff.)

Und nun der Lohn: Interpolieren ist genau dieselbe gewichtete Summe wie in 1D, nur mit drei Gewichten. Kennen wir die Temperaturen an den drei Ecken — sagen wir \(T_1 = 20\ °\mathrm{C}\), \(T_2 = 300\ °\mathrm{C}\), \(T_3 = 140\ °\mathrm{C}\) —, dann ist die Temperatur bei \(P\):

\[ T(P) = \lambda_1 T_1 + \lambda_2 T_2 + \lambda_3 T_3 = 0{,}5\cdot 20 + 0{,}25\cdot 300 + 0{,}25\cdot 140 = 120\ °\mathrm{C} . \]

Genau diese Rechnung gießen wir in wenige Zeilen reines Python. Der Code sagt Wort für Wort, was oben auf Papier stand — Fläche als halbe Kreuzdifferenz, Anteile als Verhältnis, Interpolation als gewichtete Summe:

Code

def dreiecksflaeche(p1, p2, p3):
    """Flaeche eines Dreiecks aus drei Ecken (halbe Kreuzdifferenz)."""
    kreuz = (p2[0] - p1[0]) * (p3[1] - p1[1]) \
        - (p3[0] - p1[0]) * (p2[1] - p1[1])
    if kreuz < 0.0:
        kreuz = -kreuz
    return kreuz / 2.0


def flaechenanteile(ecken, punkt):
    """Die drei Flaechenanteile eines Punktes im Dreieck (Summe 1)."""
    gesamt = dreiecksflaeche(ecken[0], ecken[1], ecken[2])
    anteil_1 = dreiecksflaeche(punkt, ecken[1], ecken[2]) / gesamt
    anteil_2 = dreiecksflaeche(punkt, ecken[2], ecken[0]) / gesamt
    anteil_3 = dreiecksflaeche(punkt, ecken[0], ecken[1]) / gesamt
    return [anteil_1, anteil_2, anteil_3]


def interpoliere(ecken, eckwerte, punkt):
    """Wert an einem Punkt aus den drei Eckwerten, gewichtet mit den Anteilen."""
    anteile = flaechenanteile(ecken, punkt)
    wert = 0.0
    for ecke in range(3):
        wert = wert + eckwerte[ecke] * anteile[ecke]
    return wert


ECKEN = [(0.0, 0.0), (4.0, 0.0), (0.0, 4.0)]
ECKWERTE = [20.0, 300.0, 140.0]
PUNKT = (1.0, 1.0)

anteile = flaechenanteile(ECKEN, PUNKT)
print("Flaechenanteile:")
for ecke in range(3):
    print("  Ecke", ecke + 1, ":", round(anteile[ecke], 3))
print("Summe der Anteile:", round(anteile[0] + anteile[1] + anteile[2], 3))
print("Interpolierte Temperatur bei P:", interpoliere(ECKEN, ECKWERTE, PUNKT),
      "Grad")

Flaechenanteile:
  Ecke 1 : 0.5
  Ecke 2 : 0.25
  Ecke 3 : 0.25
Summe der Anteile: 1.0
Interpolierte Temperatur bei P: 120.0 Grad

Interpretation der Ausgabe: Genau die Zahlen der Handrechnung — \((0{,}5;\ 0{,}25;\ 0{,}25)\), Summe 1, und 120 °C bei \(P\). Der Code hat nichts Neues getan, er hat die Papierrechnung nur nachgesprochen. Das ist der einzige neue Baustein des ganzen Kapitels; alles Weitere kennst du schon.

10.3 Drei Punkte zum Anfassen: Ecke, Kante, Mitte

Damit die Flächenanteile Anschauung bekommen, schauen wir uns drei besondere Lagen von \(P\) an. Sie zeigen, warum die Anteile genau die richtigen Gewichte für die Interpolation sind.

Code

import matplotlib.pyplot as plt

V1 = (0.0, 0.0)
V2 = (4.0, 0.0)
V3 = (0.0, 4.0)

lagen = [
    ("P auf Ecke 1", V1, "(1, 0, 0)"),
    ("P auf Kantenmitte 1–2", (2.0, 0.0), "(½, ½, 0)"),
    ("P im Schwerpunkt", (4.0 / 3.0, 4.0 / 3.0), "(⅓, ⅓, ⅓)"),
]

fig, achsen = plt.subplots(1, 3, figsize=(11.0, 3.8))
for ax, (titel, P, tripel) in zip(achsen, lagen):
    for A, B in ((V1, V2), (V2, V3), (V3, V1)):
        ax.plot([A[0], B[0]], [A[1], B[1]], color="black", lw=1.6)
    for V in (V1, V2, V3):
        ax.plot([P[0], V[0]], [P[1], V[1]], color="0.6", lw=0.8)
    ax.plot(*P, "o", color="tab:red", ms=11, zorder=5)
    ax.text(-0.3, -0.5, "1", fontsize=10, color="tab:blue")
    ax.text(4.0, -0.5, "2", fontsize=10, color="tab:blue")
    ax.text(-0.3, 4.1, "3", fontsize=10, color="tab:blue")
    ax.set_title(titel + "\n" + tripel, fontsize=10)
    ax.set_xlim(-0.8, 4.6)
    ax.set_ylim(-0.9, 4.6)
    ax.set_aspect("equal")
    ax.axis("off")
plt.tight_layout()
plt.show()

Abbildung 10.4: Drei Lagen des Punktes P und ihre Anteil-Tripel. Auf einer Ecke gehört alles dieser Ecke (1, 0, 0). Auf der Mitte einer Kante teilen sich die beiden Kantenecken hälftig, die Gegenecke zählt gar nicht (½, ½, 0). Im Schwerpunkt sind alle drei gleich (⅓, ⅓, ⅓). Das mittlere Bild ist der Schlüssel zur Kantenkompatibilität.

Was man hier sieht: drei sprechende Fälle.

Auf Ecke 1 ist der Anteil \((1, 0, 0)\): Die Interpolation gibt genau \(T_1\) zurück. Das ist die entscheidende Eigenschaft jeder Ansatzfunktion — an ihrem eigenen Knoten voll, an den anderen null. Es ist dieselbe Regel wie beim Hütchen (Kapitel 7).
Auf der Kantenmitte 1–2 ist der Anteil \((\tfrac12, \tfrac12, 0)\): Die Gegenecke 3 zählt gar nicht. Der Wert dort ist der schlichte Mittelwert von \(T_1\) und \(T_2\).
Im Schwerpunkt ist alles gleich, \((\tfrac13, \tfrac13, \tfrac13)\) — der Mittelwert aller drei Ecken.

Der mittlere Fall verdient besondere Aufmerksamkeit, denn er löst ein Problem, das man leicht übersieht: Wie passen benachbarte Dreiecke zusammen?

Kantenkompatibilität: warum das Netz nahtlos ist

Zwei Dreiecke, die sich eine Kante teilen, haben auf dieser Kante dieselben zwei Eckknoten — und die Gegenecke des jeweils anderen Dreiecks zählt auf der Kante nicht (Anteil 0). Deshalb hängt die Temperatur entlang der gemeinsamen Kante nur von den beiden geteilten Knotenwerten ab — für beide Dreiecke gleich.

Die Folge: Zwischen zwei Dreiecken gibt es keinen Sprung und keine Überlappung. Der Streckenzug von Kapitel 7 wird in der Fläche zu einer zusammenhängenden, geknickten „Zeltlandschaft” ohne Risse. Genau das brauchen wir, damit die vielen kleinen Dreiecke ein sinnvolles Ganzes ergeben — und darauf baut Kapitel 11 auf, wenn Rechtecke die Dreiecke ablösen.

10.3.1 Der Punkt wandert: die Teilflächen atmen

Wenn \(P\) durch das Dreieck wandert, ändern sich die drei Teilflächen laufend — die eine wächst, während die andere schrumpft —, ihre Summe bleibt aber immer die volle Dreiecksfläche. Und der interpolierte Wert läuft mit. Die folgende Animation zeigt das: \(P\) zieht seine Bahn, die drei Teildreiecke atmen, und die Temperatur bei \(P\) wandert zwischen den Eckwerten hin und her.

Code

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

V1 = (0.0, 0.0); V2 = (4.0, 0.0); V3 = (0.0, 4.0)
ECKWERTE = [20.0, 300.0, 140.0]
gesamt = 8.0

def flaeche(a, b, c):
    return abs((b[0]-a[0])*(c[1]-a[1]) - (c[0]-a[0])*(b[1]-a[1])) / 2.0

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10.0, 4.2),
                               gridspec_kw={"width_ratios": [1.5, 1]})

def zeichne(frame):
    ax1.clear(); ax2.clear()
    t = 2 * np.pi * frame / 36.0
    P = (4.0/3.0 + 0.9*np.cos(t), 4.0/3.0 + 0.9*np.sin(t))
    a1 = flaeche(P, V2, V3) / gesamt
    a2 = flaeche(P, V3, V1) / gesamt
    a3 = flaeche(P, V1, V2) / gesamt
    T = a1*ECKWERTE[0] + a2*ECKWERTE[1] + a3*ECKWERTE[2]
    ax1.fill([P[0], V2[0], V3[0]], [P[1], V2[1], V3[1]], color="tab:blue", alpha=0.45)
    ax1.fill([P[0], V3[0], V1[0]], [P[1], V3[1], V1[1]], color="tab:orange", alpha=0.45)
    ax1.fill([P[0], V1[0], V2[0]], [P[1], V1[1], V2[1]], color="tab:green", alpha=0.45)
    ax1.plot([V1[0],V2[0],V3[0],V1[0]], [V1[1],V2[1],V3[1],V1[1]], "k-", lw=1.8)
    ax1.plot(*P, "o", color="black", ms=9, zorder=5)
    ax1.set_xlim(-0.7, 4.5); ax1.set_ylim(-0.7, 4.5)
    ax1.set_aspect("equal"); ax1.axis("off")
    ax1.set_title("T bei P = %.0f °C" % T, fontsize=12)
    ax2.bar(["1", "2", "3"], [a1, a2, a3],
            color=["tab:blue", "tab:orange", "tab:green"])
    ax2.set_ylim(0, 1.0)
    ax2.set_ylabel("Flächenanteil")
    ax2.set_title("Anteile (Summe 1)", fontsize=10)
    return []

ani = FuncAnimation(fig, zeichne, frames=36, interval=120, blit=False)
plt.close(fig)
HTML(ani.to_jshtml())

Abbildung 10.5: Der wandernde Punkt. P zieht eine Schleife durchs Dreieck; die drei Teilflächen atmen (mal groß, mal klein), und die drei Anteile — als Balken rechts — verschieben sich entsprechend. Der interpolierte Temperaturwert oben läuft zwischen den Eckwerten 20, 300 und 140 °C mit. Nirgends wird ein Anteil negativ, weil P im Dreieck bleibt.

Was man sieht: Die Teilflächen sind kommunizierende Röhren — was die eine gewinnt, verliert eine andere, die Summe bleibt konstant. Deshalb summieren sich die Anteile immer zu 1, und der interpolierte Wert bleibt stets zwischen dem kleinsten und größten Eckwert. Springt \(P\) nahe an eine Ecke, dominiert deren Anteil und die Temperatur nähert sich deren Eckwert.

10.4 Das Zelt: die Ansatzfunktion in der Fläche

Jetzt der Brückenschlag zu Kapitel 7. Dort war die Ansatzfunktion eines Knotens ein Hütchen: über dem eigenen Knoten auf 1, zu den Nachbarn geradlinig auf 0. In der Fläche wird daraus ein Zelt: Betrachte den Anteil \(\lambda_k\) eines Knotens als Funktion des Ortes — er ist 1 am eigenen Knoten und fällt über jedes anliegende Dreieck geradlinig auf 0 an den anderen Ecken. Setzt man die schrägen Dreiecksflächen aller Dreiecke zusammen, die an diesem Knoten hängen, entsteht ein Zeltdach.

Code

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri

# Sechs Dreiecke fächern um einen Mittelknoten
knoten = [(0.0, 0.0)]
for k in range(6):
    winkel = 2 * np.pi * k / 6.0
    knoten.append((np.cos(winkel), np.sin(winkel)))
dreiecke = []
for k in range(6):
    dreiecke.append((0, 1 + k, 1 + (k + 1) % 6))
xs = [p[0] for p in knoten]
ys = [p[1] for p in knoten]
# Zelt-Hoehe: 1 am Mittelknoten, 0 am Rand
zs = [1.0] + [0.0] * 6
triang = mtri.Triangulation(xs, ys, dreiecke)

fig = plt.figure(figsize=(11.0, 4.2))
ax1 = fig.add_subplot(1, 2, 1, projection="3d")
ax1.plot_trisurf(triang, zs, cmap="viridis", edgecolor="0.3", linewidth=0.4)
ax1.set_title("Das Zelt (3D)", fontsize=10)
ax1.set_zlim(0, 1.1)
ax1.set_xlabel("x"); ax1.set_ylabel("y")

ax2 = fig.add_subplot(1, 2, 2)
ax2.tricontourf(triang, zs, levels=10, cmap="viridis")
hoehen = ax2.tricontour(triang, zs, levels=6, colors="white", linewidths=0.8)
ax2.triplot(triang, color="0.7", lw=0.6)
ax2.plot(0, 0, "o", color="red", ms=9)
ax2.set_title("Draufsicht mit Höhenlinien", fontsize=10)
ax2.set_aspect("equal")
ax2.set_xlabel("x"); ax2.set_ylabel("y")
plt.tight_layout()
plt.show()

Abbildung 10.6: Das Zelt eines Knotens über seinen sechs Dreiecken — die flächenhafte Fortsetzung des Hütchens aus Kapitel 7. Links die 3D-Ansicht: über dem eigenen Knoten steht das Zelt auf der Höhe 1, an allen Nachbarknoten auf 0, dazwischen jede Dreiecksfläche eben (schräg). Rechts die Draufsicht mit Höhenlinien: sechseckige Ringe, die zum Zentrum hin auf 1 steigen.

Was man hier sieht: links das Zelt in 3D — ein spitzes Dach, das über dem Mittelknoten auf 1 steht und zu allen sechs Nachbarn auf 0 abfällt; jede einzelne Dreiecksfläche ist eben (das garantieren die drei Ecken, die eine Ebene festlegen). Rechts dieselbe Funktion von oben, mit Höhenlinien: konzentrische sechseckige Ringe, die zum roten Mittelknoten hin auf 1 steigen. Genau wie in 1D überlappen sich nur die Zelte benachbarter Knoten, und an jedem Punkt der Fläche summieren sich alle Zelte zu 1.

Das Hütchen, eine Dimension höher

Schneidet man das Zelt mit einer senkrechten Ebene durch den Mittelknoten, erhält man exakt das Hütchen aus Kapitel 7: 1 in der Mitte, geradlinig auf 0 zu beiden Seiten. Das Zelt ist das Hütchen, in die zweite Dimension gezogen. Und weil jede Zeltfläche eben ist, hat sie über jedem Dreieck eine konstante Steigung — in \(x\)-Richtung eine, in \(y\)-Richtung eine. Das ist das flächenhafte Gegenstück zu den \(\pm 1/h\) des Hütchens, und es macht die Elementmatrix gleich genauso einfach wie in 1D.

10.5 Das Netz als Bauplan

Bevor wir rechnen, brauchen wir eine Beschreibung des Netzes, die ein Programm lesen kann. Sie besteht aus zwei Listen — genau den Datenstrukturen aus Kapitel 2:

die Knotenliste: für jeden Knoten seine Koordinaten \((x, y)\),
die Dreiecksliste: für jedes Dreieck die drei Knotennummern seiner Ecken.

Mehr ist ein Netz nicht. Man liest es wie einen Bauplan: Die Dreiecksliste sagt, welche Knoten ein Dreieck bilden, die Knotenliste sagt, wo diese Knoten liegen. Abbildung 10.7 zeigt beide Listen neben dem gezeichneten Netz.

Code

import matplotlib.pyplot as plt

knoten = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0),
          (0.0, 0.05), (0.5, 0.05), (1.0, 0.05),
          (0.0, 0.10), (0.5, 0.10), (1.0, 0.10)]
dreiecke = [(0, 1, 4), (0, 4, 3), (1, 2, 5), (1, 5, 4),
            (3, 4, 7), (3, 7, 6), (4, 5, 8), (4, 8, 7)]

fig, (axl, axr) = plt.subplots(1, 2, figsize=(11.5, 4.2),
                               gridspec_kw={"width_ratios": [1, 1.4]})

# Linke Seite: die zwei Listen als Text
zeilen = ["Knotenliste (Nr: x, y):"]
for i, (x, y) in enumerate(knoten):
    zeilen.append("  %d: (%.2f, %.2f)" % (i + 1, x, y))
zeilen.append("")
zeilen.append("Dreiecksliste (Ecken):")
for t, tri in enumerate(dreiecke):
    zeilen.append("  T%d: %d – %d – %d" % (t + 1, tri[0]+1, tri[1]+1, tri[2]+1))
axl.text(0.0, 1.0, "\n".join(zeilen), family="monospace", fontsize=9,
         va="top")
axl.axis("off")

# Rechte Seite: das gezeichnete Netz
bx = [p[0] for p in knoten]
by = [p[1] for p in knoten]
axr.triplot(bx, by, dreiecke, color="0.6", lw=1.2)
axr.plot(bx, by, "o", color="tab:orange", ms=16, markeredgecolor="black",
         zorder=4)
for i, (x, y) in enumerate(knoten):
    axr.text(x, y, "%d" % (i + 1), ha="center", va="center", fontsize=8,
             zorder=5, weight="bold")
# Dreiecksnummern in den Schwerpunkten
for t, tri in enumerate(dreiecke):
    sx = (knoten[tri[0]][0] + knoten[tri[1]][0] + knoten[tri[2]][0]) / 3
    sy = (knoten[tri[0]][1] + knoten[tri[1]][1] + knoten[tri[2]][1]) / 3
    axr.text(sx, sy, "T%d" % (t + 1), ha="center", va="center", fontsize=8,
             color="tab:blue", style="italic")
axr.set_xlabel("Länge x (m)")
axr.set_ylabel("Höhe y (m)")
axr.set_title("9 Knoten, 8 Dreiecke", fontsize=10)
axr.set_aspect("equal")
plt.tight_layout()
plt.show()

Abbildung 10.7: Das Netz des Trägers als Bauplan. Links die beiden Listen — neun Knoten mit Koordinaten, acht Dreiecke mit je drei Knotennummern. Rechts das gezeichnete Netz mit nummerierten Knoten (orange) und Dreiecken (grau, eingekreist). Die Knoten sind von unten nach oben, die Dreiecke zellenweise durchnummeriert; jede Gitterzelle ist von unten links nach oben rechts geteilt.

Was man hier sieht: links der komplette Bauplan in zwei Listen, rechts das Netz. Die neun Knoten sitzen auf einem \(3\times 3\)-Gitter (\(x = 0 / 0{,}5 / 1\ \mathrm{m}\), \(y = 0 / 0{,}05 / 0{,}1\ \mathrm{m}\)), von unten nach oben nummeriert. Knoten 5 sitzt genau in der Mitte des Trägers. Jede der vier Gitterzellen ist diagonal in zwei Dreiecke geteilt — macht acht. Diese beiden Listen sind alles, was das Programm über die Geometrie wissen muss.

10.6 Die Elementmatrix: der 3×3-Steckbrief

Jetzt der Schritt, der in 1D die 2×2-Federmatrix lieferte — in 2D wird daraus eine 3×3-Matrix, weil ein Dreieck drei Knoten hat. Und das Rezept ist dasselbe wie in Kapitel 7: Jeder Eintrag ist Steigung mal Steigung, mal Material, mal Fläche. Nur gibt es jetzt zwei Steigungen pro Zeltfläche (eine in \(x\), eine in \(y\)), und wir zählen ihre Produkte zusammen:

\[ K_{ij} = \lambda \cdot d \cdot A \cdot \Bigl( \underbrace{\tfrac{\partial N_i}{\partial x}\,\tfrac{\partial N_j}{\partial x}}_{\text{Steigungen in }x} \;+\; \underbrace{\tfrac{\partial N_i}{\partial y}\,\tfrac{\partial N_j}{\partial y}}_{\text{Steigungen in }y} \Bigr) . \]

In Worten: Leitfähigkeit \(\lambda\) mal Tiefe \(d\) mal Dreiecksfläche \(A\) mal (Steigung mal Steigung, in beiden Richtungen zusammengezählt). Weil die Steigungen über das ganze Dreieck konstant sind, ist der Klammerausdruck konstant — und das Elementintegral wieder nur „Wert mal Fläche”, genau wie in 1D. Die Steigungen der Zeltfläche von Knoten \(i\) bekommt man aus den Ecken:

\[ \frac{\partial N_i}{\partial x} = \frac{y_j - y_k}{2A}, \qquad \frac{\partial N_i}{\partial y} = \frac{x_k - x_j}{2A}, \]

wobei \(j, k\) die beiden anderen Ecken sind und \(2A\) die doppelte Fläche. Das sieht abstrakt aus — machen wir es an einem echten Dreieck des Trägers konkret.

10.6.1 Ein Dreieck von Hand

Wir nehmen das rechtwinklige Dreieck T1 aus dem Netz: die Knoten 1, 2, 5 mit den Ecken \((0,0)\), \((0{,}5, 0)\), \((0{,}5, 0{,}05)\). Rechtwinklig heißt: Die Zahlen werden freundlich. Zuerst die doppelte Fläche und die Fläche selbst:

\[ 2A = 0{,}025\ \mathrm{m^2}, \qquad A = 0{,}0125\ \mathrm{m^2} \]

(ein rechtwinkliges Dreieck mit Katheten \(0{,}5\) und \(0{,}05\): Fläche \(\tfrac12 \cdot 0{,}5 \cdot 0{,}05 = 0{,}0125\)). Nun die Steigungen der drei Zeltflächen — je eine in \(x\), eine in \(y\), in \(1/\mathrm{m}\):

Knoten	\(\partial N/\partial x\)	\(\partial N/\partial y\)
1 (0, 0)	\(-2\)	\(0\)
2 (0,5, 0)	\(+2\)	\(-20\)
5 (0,5, 0,05)	\(0\)	\(+20\)

Die Steigungen summieren sich spaltenweise zu null (das müssen sie, weil die drei Zelte zusammen überall die konstante Höhe 1 ergeben). Die \(y\)-Steigungen sind mit \(\pm 20\) groß, weil das Dreieck nur \(0{,}05\ \mathrm{m}\) hoch ist — auf kurzer Höhe muss das Zelt steil abfallen. Der gemeinsame Vorfaktor ist

\[ \lambda \cdot d \cdot A = 50 \cdot 0{,}1 \cdot 0{,}0125 = 0{,}0625\ \tfrac{\mathrm{W}}{\mathrm{K}} . \]

Jetzt jeder Eintrag als „Steigung mal Steigung, mal \(0{,}0625\)“. Zwei Beispiele, den Rest macht das Programm:

\[ \begin{aligned} K_{11} &= 0{,}0625 \cdot \bigl((-2)(-2) + (0)(0)\bigr) = 0{,}0625 \cdot 4 = 0{,}25,\\ K_{25} &= 0{,}0625 \cdot \bigl((+2)(0) + (-20)(+20)\bigr) = 0{,}0625 \cdot (-400) = -25 . \end{aligned} \]

Alle neun zusammengeschrieben ergeben den Steckbrief von T1 (in W/K):

\[ \mathbf{K}_{\text{T1}} = \begin{pmatrix} 0{,}25 & -0{,}25 & 0\\ -0{,}25 & 25{,}25 & -25\\ 0 & -25 & 25 \end{pmatrix} \qquad \substack{\text{Zeilen 1, 2, 5}} \]

Und hier die vertraute Probe: Jede Zeile summiert sich zu null (\(0{,}25-0{,}25+0=0\), \(-0{,}25+25{,}25-25=0\), \(0-25+25=0\)). Das ist dasselbe Zeichen wie bei der Feder in Kapitel 3 — verschiebt man alle drei Knoten um dasselbe Stück (oder erhitzt sie gleich), passiert nichts, also müssen die Zeilen sich wegheben. Lassen wir das Programm dieselbe Matrix ausrechnen, damit die Handrechnung belegt ist:

Code

import os
import sys
sys.path.insert(0, os.path.join("..", "programme", "gemeinsam"))
from dreiecke import baue_dreieck_matrix

ecken_t1 = [(0.0, 0.0), (0.5, 0.0), (0.5, 0.05)]
steckbrief = baue_dreieck_matrix(ecken_t1, 50.0, 0.1)

print("Steckbrief von T1 (Knoten 1, 2, 5), in W/K:")
for zeile in steckbrief:
    zeilensumme = zeile[0] + zeile[1] + zeile[2]
    print("  [ %7.2f %7.2f %7.2f ]   Zeilensumme = %.2f"
          % (zeile[0], zeile[1], zeile[2], zeilensumme))

Steckbrief von T1 (Knoten 1, 2, 5), in W/K:
  [    0.25   -0.25    0.00 ]   Zeilensumme = 0.00
  [   -0.25   25.25  -25.00 ]   Zeilensumme = 0.00
  [    0.00  -25.00   25.00 ]   Zeilensumme = 0.00

Interpretation der Ausgabe: Zahl für Zahl die Matrix der Handrechnung, und jede Zeilensumme ist null. Die Handrechnung ist damit belegt. Ein neuer Elementtyp, dieselbe Prüfregel wie seit Kapitel 3.

10.7 Assemblieren, einspannen, lösen

Ab hier ist es reine Wiederholung — das vierte Mal in diesem Buch, dass wir Element-Steckbriefe zu einem System zusammenschieben. Jedes der acht Dreiecke legt seinen 3×3-Steckbrief an die Stelle seiner drei Knoten in die große \(9\times 9\)-Matrix, und wo Dreiecke sich einen Knoten teilen, addieren sich ihre Beiträge auf dessen Diagonalfeld. Die folgende Animation pflastert den Träger Dreieck für Dreieck und lässt das System mitwachsen.

Code

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

knoten = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0),
          (0.0, 0.05), (0.5, 0.05), (1.0, 0.05),
          (0.0, 0.10), (0.5, 0.10), (1.0, 0.10)]
dreiecke = [(0, 1, 4), (0, 4, 3), (1, 2, 5), (1, 5, 4),
            (3, 4, 7), (3, 7, 6), (4, 5, 8), (4, 8, 7)]

def elementmatrix(tri):
    (xa,ya),(xb,yb),(xc,yc) = knoten[tri[0]],knoten[tri[1]],knoten[tri[2]]
    zwei = xa*(yb-yc)+xb*(yc-ya)+xc*(ya-yb)
    fl = abs(zwei)/2.0
    sx = [(yb-yc)/zwei,(yc-ya)/zwei,(ya-yb)/zwei]
    sy = [(xc-xb)/zwei,(xa-xc)/zwei,(xb-xa)/zwei]
    K = [[50.0*0.1*fl*(sx[i]*sx[j]+sy[i]*sy[j]) for j in range(3)] for i in range(3)]
    return K

def matrix_nach(anzahl):
    M = [[0.0]*9 for _ in range(9)]
    for e in range(anzahl):
        tri = dreiecke[e]; K = elementmatrix(tri)
        for a in range(3):
            for b in range(3):
                M[tri[a]][tri[b]] += K[a][b]
    return M

bx = [p[0] for p in knoten]; by = [p[1] for p in knoten]
fig, (axl, axr) = plt.subplots(1, 2, figsize=(11.0, 4.4),
                               gridspec_kw={"width_ratios": [1.3, 1]})

def zeichne(schritt):
    axl.clear(); axr.clear()
    # linkes Netz
    for e in range(schritt):
        tri = dreiecke[e]
        xs = [knoten[tri[0]][0], knoten[tri[1]][0], knoten[tri[2]][0]]
        ys = [knoten[tri[0]][1], knoten[tri[1]][1], knoten[tri[2]][1]]
        farbe = "tab:orange" if e == schritt-1 else "#cfe3f5"
        axl.fill(xs, ys, color=farbe, edgecolor="0.4")
    axl.plot(bx, by, "o", color="tab:orange", ms=9, markeredgecolor="black")
    axl.set_xlim(-0.1, 1.1); axl.set_ylim(-0.03, 0.13)
    axl.set_aspect("equal"); axl.axis("off")
    axl.set_title("Dreiecke gepflastert: %d von 8" % schritt, fontsize=10)
    # rechte Matrix
    M = matrix_nach(schritt)
    aktiv = dreiecke[schritt-1] if schritt >= 1 else ()
    for r in range(9):
        for c in range(9):
            gehoert = r in aktiv and c in aktiv
            fl = "#ffe0b3" if gehoert else "white"
            axr.add_patch(plt.Rectangle((c, 8-r), 1, 1, facecolor=fl,
                                        edgecolor="0.8", lw=0.6))
            if abs(M[r][c]) > 1e-9:
                axr.add_patch(plt.Rectangle((c, 8-r), 1, 1, facecolor="tab:blue",
                                            alpha=0.35, edgecolor="0.8", lw=0.6))
    axr.set_xlim(0, 9); axr.set_ylim(0, 9)
    axr.set_aspect("equal"); axr.axis("off")
    axr.set_title("9×9-System (belegt = blau)", fontsize=10)
    return []

ani = FuncAnimation(fig, zeichne, frames=9, interval=900, blit=False)
plt.close(fig)
HTML(ani.to_jshtml())

Abbildung 10.8: Die Assemblierung des 2D-Systems — die vierte Wiederkehr des Schiebespiels aus Kapitel 3. Links wird der Träger Dreieck für Dreieck gepflastert (das jeweils neue Dreieck orange), rechts füllt sich die 9×9-Systemmatrix mit. Wo ein Knoten zu mehreren Dreiecken gehört, wächst seine Diagonalzahl mit jedem hinzukommenden Dreieck. Am Ende steht das ganze System.

Was man sieht: Der Träger füllt sich mit Dreiecken, und im gleichen Takt füllt sich die Matrix. Anders als in 1D, wo jeder Knoten nur zwei Nachbarn hatte (Tridiagonale), hängt ein Knoten hier mit mehreren zusammen — auch schräg über die Netzkanten. Das Band ist deshalb breiter als in 1D, bleibt aber ein Band: weit weg von der Diagonale ist alles leer.

10.7.1 Das System, ausgeschrieben

Schreiben wir die assemblierte Matrix als Textblock — die Matrix-Drucker-Konvention aus Kapitel 3, Nullen als Punkte. Daneben zum Vergleich das tridiagonale 1D-System aus Kapitel 4/Kapitel 7.

Code

import os
import sys
sys.path.insert(0, os.path.join("..", "programme", "gemeinsam"))
from dreiecke import assembliere_2d
from matrix_drucker import drucke_matrix

knoten = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0),
          (0.0, 0.05), (0.5, 0.05), (1.0, 0.05),
          (0.0, 0.10), (0.5, 0.10), (1.0, 0.10)]
dreiecke = [(0, 1, 4), (0, 4, 3), (1, 2, 5), (1, 5, 4),
            (3, 4, 7), (3, 7, 6), (4, 5, 8), (4, 8, 7)]

matrix = assembliere_2d(knoten, dreiecke, 50.0, 0.1)
print("Systemmatrix des 8-Dreiecke-Trägers (9x9, W/K), Nullen als Punkte:")
print()
drucke_matrix(matrix)

Systemmatrix des 8-Dreiecke-Trägers (9x9, W/K), Nullen als Punkte:

      25      -0       .     -25       .       .       .       .       .
      -0      50      -0       .     -50       .       .       .       .
       .      -0      25       .       .     -25       .       .       .
     -25       .       .      50      -0       .     -25       .       .
       .     -50       .      -0     101      -0       .     -50       .
       .       .     -25       .      -0      50       .       .     -25
       .       .       .     -25       .       .      25      -0       .
       .       .       .       .     -50       .      -0      50      -0
       .       .       .       .       .     -25       .      -0      25

Interpretation der Ausgabe: ein breiteres Band als die 1D-Tridiagonale. Der Matrix-Drucker rundet auf ganze W/K. Auf der Diagonale stehen die großen Zahlen (gerundet 25, 50 und 101; in Wahrheit 25,25 und 50,5), daneben die starke Kopplung quer über die kurze Höhe (\(-25\) und \(-50\)). Die Kopplung längs des Trägers ist dagegen so schwach — nur \(-0{,}25\) bzw. \(-0{,}5\ \mathrm{W/K}\) —, dass der ganzzahlige Drucker sie als „\(-0\)” zeigt: winzig, aber eben nicht null (einen Punkt druckt er nur, wo gar keine Kopplung ist). Die lang gestreckten, flachen Dreiecke koppeln über die kurze Höhe viel stärker als über die lange Länge — die Geometrie steht der Matrix ins Gesicht geschrieben. Abbildung 10.9 stellt beide Bandstrukturen nebeneinander.

Code

import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11.0, 4.4))

# 1D: 5x5 tridiagonal
n1 = 5
for r in range(n1):
    for c in range(n1):
        belegt = abs(r - c) <= 1
        farbe = "tab:orange" if belegt else "white"
        ax1.add_patch(plt.Rectangle((c, n1-1-r), 1, 1, facecolor=farbe,
                                    edgecolor="0.7"))
ax1.set_xlim(0, n1); ax1.set_ylim(0, n1)
ax1.set_aspect("equal"); ax1.axis("off")
ax1.set_title("1D-Träger (5×5): tridiagonal", fontsize=10)

# 2D: 9x9 aus der echten Matrix
knoten = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0),
          (0.0, 0.05), (0.5, 0.05), (1.0, 0.05),
          (0.0, 0.10), (0.5, 0.10), (1.0, 0.10)]
dreiecke = [(0, 1, 4), (0, 4, 3), (1, 2, 5), (1, 5, 4),
            (3, 4, 7), (3, 7, 6), (4, 5, 8), (4, 8, 7)]
import sys, os
sys.path.insert(0, os.path.join("..", "programme", "gemeinsam"))
from dreiecke import assembliere_2d
M = assembliere_2d(knoten, dreiecke, 50.0, 0.1)
n2 = 9
for r in range(n2):
    for c in range(n2):
        belegt = abs(M[r][c]) > 1e-9
        farbe = "tab:blue" if belegt else "white"
        ax2.add_patch(plt.Rectangle((c, n2-1-r), 1, 1, facecolor=farbe,
                                    edgecolor="0.8"))
ax2.set_xlim(0, n2); ax2.set_ylim(0, n2)
ax2.set_aspect("equal"); ax2.axis("off")
ax2.set_title("2D-Träger (9×9): breiteres Band", fontsize=10)
plt.tight_layout()
plt.show()

Abbildung 10.9: Zwei Bandstrukturen. Links das 1D-System aus Kapitel 7 (5×5, tridiagonal): nur die Diagonale und je ein Nachbar. Rechts das 2D-System dieses Kapitels (9×9): dasselbe Prinzip — belegt nur nahe der Diagonale —, aber ein breiteres Band, weil jeder Knoten in der Fläche mehr Nachbarn hat. Beide sind dünn besetzt, beide löst derselbe Löser aus Kapitel 4.

Was man hier sieht: links das schmale 1D-Band, rechts das breitere 2D-Band — aber beide sind dünn besetzt, beide haben ihre Zahlen dicht um die Diagonale. Der Unterschied ist nur die Bandbreite, nicht das Prinzip. Und weil beide dünn besetzt sind, löst sie derselbe Löser aus Kapitel 4.

10.7.2 Randwerte und Lösung

Die linke Kante (Knoten 1, 4, 7) halten wir auf 20 °C, die rechte (Knoten 3, 6, 9) auf 300 °C — das Streich-Rezept aus Kapitel 8, jetzt auf sechs statt zwei Randknoten angewandt. Übrig bleiben die drei freien Knoten der Mittelspalte (2, 5, 8). Gelöst wird mit Gauß-Seidel aus Kapitel 4.

Vorhersage-Punkt

Bevor wir lösen: Knoten 5 sitzt genau in der Trägermitte. Links wird der Träger auf 20 °C gehalten, rechts auf 300 °C. Welche Temperatur erwartest du an Knoten 5? Und was ist mit den anderen beiden Mittelknoten (2 unten, 8 oben) — gleiche Temperatur oder anders? Leg dich fest.

Code

import os
import sys
sys.path.insert(0, os.path.join("..", "programme", "gemeinsam"))
from dreiecke import assembliere_2d, reduziere_system, setze_ein
from loeser import gauss_seidel

knoten = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0),
          (0.0, 0.05), (0.5, 0.05), (1.0, 0.05),
          (0.0, 0.10), (0.5, 0.10), (1.0, 0.10)]
dreiecke = [(0, 1, 4), (0, 4, 3), (1, 2, 5), (1, 5, 4),
            (3, 4, 7), (3, 7, 6), (4, 5, 8), (4, 8, 7)]

matrix = assembliere_2d(knoten, dreiecke, 50.0, 0.1)

# Randwerte: linke Kante 20, rechte 300, Mittelspalte frei (None).
randwerte = [20.0, None, 300.0, 20.0, None, 300.0, 20.0, None, 300.0]
kleine_matrix, kleine_seite, frei = reduziere_system(
    matrix, [0.0] * 9, randwerte)

startwerte = [0.0] * len(frei)
loesung_frei = gauss_seidel(kleine_matrix, kleine_seite, startwerte,
                            0.00000001, 2000)
temperaturen = setze_ein(loesung_frei, frei, randwerte)

frei_namen = []
for k in frei:
    frei_namen.append(k + 1)
print("Freie Knoten (Mittelspalte):", frei_namen)
print()
print("Knoten |   Ort (x, y)    | Temperatur")
print("-------+-----------------+-----------")
for i in range(9):
    x, y = knoten[i]
    print("   %d   | (%.2f, %.3f) m | %6.1f Grad" % (i + 1, x, y, temperaturen[i]))

Freie Knoten (Mittelspalte): [2, 5, 8]

Knoten |   Ort (x, y)    | Temperatur
-------+-----------------+-----------
   1   | (0.00, 0.000) m |   20.0 Grad
   2   | (0.50, 0.000) m |  160.0 Grad
   3   | (1.00, 0.000) m |  300.0 Grad
   4   | (0.00, 0.050) m |   20.0 Grad
   5   | (0.50, 0.050) m |  160.0 Grad
   6   | (1.00, 0.050) m |  300.0 Grad
   7   | (0.00, 0.100) m |   20.0 Grad
   8   | (0.50, 0.100) m |  160.0 Grad
   9   | (1.00, 0.100) m |  300.0 Grad

Interpretation der Ausgabe: Alle drei Mittelknoten stehen auf 160 °C — genau in der Mitte zwischen 20 und 300, und oben wie unten gleich. Lag deine Vorhersage richtig? Das ist die Pointe des Kapitels:

Vermutung aufgelöst

Die Vermutung war: „2D braucht ganz neue Mathematik.” Und heraus kommt: dieselbe Gerade wie in 1D. Bei diesem Rand (links überall kühl am Auflager, rechts überall heiß beim Brandherd) hängt die Temperatur gar nicht von der Höhe \(y\) ab — sie wächst nur mit der Länge \(x\), exakt wie beim 1D-Träger aus Kapitel 7. Die Mittelspalte liegt bei \(x = 0{,}5\), also bei 160 °C. Hier ist 2D noch 1D in Verkleidung. Genau das haben wir im Kasten „Warum durften wir 1D rechnen?” vorausgesagt: Solange die Höhe gleichmäßig belastet wird, verrät sie nichts Neues. Das ändert sich, sobald die Wärme nicht mehr von der Seite, sondern von unten kommt: Treffen die Brandgase die Unterkante und hält die Betondecke die Oberseite kühl, dann kippt das Gefälle in die Höhe — ein Höhenunterschied, der kein Streifenbild mehr ist und der den Träger später (Kapitel 14, Kapitel 15) krümmt. Genau das rechnest du in Übung 10.2 nach.

Das Temperatur-Farbbild macht die Streifen sichtbar: waagerechte Bänder, kein Unterschied zwischen oben und unten.

Code

import matplotlib.pyplot as plt
import matplotlib.tri as mtri

knoten = [(0.0, 0.0), (0.5, 0.0), (1.0, 0.0),
          (0.0, 0.05), (0.5, 0.05), (1.0, 0.05),
          (0.0, 0.10), (0.5, 0.10), (1.0, 0.10)]
dreiecke = [(0, 1, 4), (0, 4, 3), (1, 2, 5), (1, 5, 4),
            (3, 4, 7), (3, 7, 6), (4, 5, 8), (4, 8, 7)]
temperaturen = [20.0, 160.0, 300.0, 20.0, 160.0, 300.0, 20.0, 160.0, 300.0]

xs = [p[0] for p in knoten]; ys = [p[1] for p in knoten]
triang = mtri.Triangulation(xs, ys, dreiecke)

fig, ax = plt.subplots(figsize=(9.0, 3.0))
bild = ax.tricontourf(triang, temperaturen, levels=20, cmap="inferno")
ax.triplot(triang, color="white", lw=0.6, alpha=0.6)
ax.plot(xs, ys, "o", color="white", ms=6, markeredgecolor="black")
for i in range(9):
    ax.annotate("%.0f" % temperaturen[i], xy=(knoten[i][0], knoten[i][1]),
                xytext=(0, 7), textcoords="offset points", ha="center",
                fontsize=8, color="white")
fig.colorbar(bild, ax=ax, label="Temperatur (°C)")
ax.set_xlabel("Länge x (m)"); ax.set_ylabel("Höhe y (m)")
ax.set_title("Temperaturfeld: waagerechte Streifen")
ax.set_aspect("equal")
plt.tight_layout()
plt.show()

Abbildung 10.10: Das Temperaturfeld des Trägers als Farbbild. Die Farbe wächst gleichmäßig von 20 °C links (dunkel) auf 300 °C rechts (hell); die Höhenlinien sind senkrecht — bei gleichem Rand hängt die Temperatur nur von der Länge ab, nicht von der Höhe. Das ist das bunte Bild zu den 160 °C der Mittelspalte.

Was man hier sieht: senkrechte Farbübergänge, waagerechte Streifen. Die Temperatur klettert gleichmäßig von links (20 °C) nach rechts (300 °C), und auf jeder senkrechten Linie ist sie konstant. Das bestätigt die 160 °C und die Pointe: Bei symmetrischem Rand ist das Flächenfeld nur die in die Breite gezogene 1D-Gerade. In Übung 10.2 kippen die Streifen: Trifft das Brandgas die Unterkante, läuft das Gefälle über die Höhe.

10.8 Die interaktive Einheit: das Dreiecks-Labor

Jetzt bist du dran. Das Dreiecks-Labor hat zwei Teile. Teil 1 lässt dich einen Punkt in ein Dreieck setzen und liest die drei Flächenanteile samt interpoliertem Wert ab. Teil 2 ist der ganze 8-Dreiecke-Träger: Randwerte setzen, System drucken, mit Gauß-Seidel lösen, Farbbild zeichnen.

10.8.1 Teil 1 — Flächenanteile eines Punktes

Setze PUNKT an eine Stelle im Dreieck und sieh, wie sich die drei Anteile und der interpolierte Wert ändern. Probiere eine Ecke (\((0,0)\)), eine Kantenmitte (\((2,0)\)) und den Schwerpunkt (etwa \((1{,}33,\ 1{,}33)\)).

Abbildung 10.11: Vorgerenderte Fassung von Teil 1 des Dreiecks-Labors: Punkt P = (1, 1) mit den Anteilen (½, ¼, ¼) und T = 120 °C. Im Browser ersetzt dein eigenes Ergebnis dieses Bild, sobald du den Punkt verschiebst.

10.8.2 Teil 2 — der 8-Dreiecke-Träger

Setze die beiden Randtemperaturen LINKS und RECHTS, und das Labor baut das System, druckt es, löst mit Gauß-Seidel und zeichnet das Farbbild. Der Code ist Zeile für Zeile das, was oben im Kapitel stand.

Abbildung 10.12: Vorgerenderte Fassung von Teil 2: das Temperaturfeld bei 20 °C links und 300 °C rechts — waagerechte Streifen, Mittelspalte 160 °C. Im Browser ersetzt dein eigenes Ergebnis dieses Bild, sobald du die Randwerte änderst.

Was man sieht: Bei gleichem Rand links und rechts bleibt das Feld streifig (Mittelspalte 160 °C). Setze in Teil 2 einmal LINKS = 300 und RECHTS = 20 — die Streifen drehen sich um. Erst Übung 10.2 (Brandgas an der Unterkante) lässt das Gefälle über die Höhe kippen und bringt so die zweite Dimension ins Spiel.

10.9 Das Kapitel-Programm

Das vollständige, eigenständig lauffähige Skript liegt in programme/kap10/kap10_dreiecke_fem.py. Es baut das 8-Dreiecke-Netz, rechnet den Steckbrief von T1 vor, druckt die assemblierte Systemmatrix, baut die Randwerte ein und löst mit Gauß-Seidel — dieselben Schritte wie im Kapitel, an einem Stück. Die 2D-Bausteine selbst (baue_dreieck_matrix, assembliere_2d, reduziere_system, …) liegen im gemeinsamen Modul programme/gemeinsam/dreiecke.py, das auch Kapitel 12 und die Verformung in Kapitel 14 wiederverwenden. Führe das Skript mit python kap10_dreiecke_fem.py aus.

10.10 Die Namensschilder

Wie in Kapitel 7 haben wir die Sache zuerst gebaut und die Fachwörter weggelassen. Jetzt die Schilder.

Die Fachwörter, jetzt nachgereicht

Baryzentrische Koordinaten: das Fachwort für unsere Flächenanteile \(\lambda_1, \lambda_2, \lambda_3\). „Baryzentrum” ist griechisch-lateinisch für Schwerpunkt (barys = schwer). Die Idee stammt von August Ferdinand Möbius (1827): Man beschreibt einen Punkt durch die Gewichte, die man auf die drei Ecken legen müsste, damit ihr Schwerpunkt genau dort liegt. Diese Gewichte sind die Flächenanteile.
P1-Dreieck (lineares Dreieckselement): das Dreieck mit einem Knoten je Ecke und linearer Interpolation dazwischen — das „P1” steht für Polynome vom Grad 1 (eben, keine Krümmung). Es ist das flächenhafte Gegenstück zum linearen Stab-Element aus Kapitel 7.
Zeltfunktion (Ansatz-/Formfunktion in 2D): die Funktion \(N_k = \lambda_k\), über dem eigenen Knoten 1, an allen anderen 0 — das Zelt, das wir gezeichnet haben. Es ist das Hütchen aus Kapitel 7, eine Dimension höher.
Triangulierung: die Zerlegung einer Fläche in Dreiecke — der Bauplan aus Knoten- und Dreiecksliste. Das zugehörige Netz ist das flächenhafte Gegenstück zur Knotenkette in 1D.

Alle diese Wörter stehen im Glossar; hier haben wir zuerst die Sache gebaut und dann das Schild geklebt.

Woher kommen die Dreiecke? Ein Namedrop.

Wir haben das Netz von Hand hingeschrieben — neun Knoten, acht Dreiecke. Für echte, krumme Formen übernimmt das ein Netzgenerator. Das bekannteste Verfahren ist die Delaunay-Triangulierung (nach Boris Delaunay, 1934): Sie verbindet die Punkte so, dass die Dreiecke möglichst „gutmütig” werden — keine extrem spitzen Nadeln, sondern nach Möglichkeit gedrungene Dreiecke. Wie das genau geht, ist ein Kapitel für sich und nicht unser Thema; wir schreiben unsere kleinen Netze weiter von Hand. Merke nur: Ein Programm kann die Dreiecksliste aus einer Punktwolke selbst erzeugen, und Delaunay ist der Klassiker dafür.

Merkkasten

Der einzige neue Baustein in der Fläche ist die Interpolation mit Flächenanteilen: Punkt \(P\) mit den Ecken verbinden, die drei Teildreiecke messen (halbe Kreuzdifferenz), durch die Gesamtfläche teilen. Die Anteile summieren sich zu 1.
Interpolieren heißt: jeder Eckwert mal seinem Anteil, aufsummiert — dieselbe gewichtete Summe wie in 1D, nur mit drei Gewichten.
Die Ansatzfunktion eines Knotens ist ein Zelt über seinen Dreiecken — das Hütchen aus Kapitel 7, eine Dimension höher. Jede Zeltfläche ist eben, also hat sie über jedem Dreieck konstante Steigungen (eine in \(x\), eine in \(y\)).
Ein Netz ist eine Knotenliste plus eine Dreiecksliste — ein Bauplan.
Die Elementmatrix ist jetzt \(3\times 3\): \(K_{ij} = \lambda\,d\,A\, (\partial_x N_i\,\partial_x N_j + \partial_y N_i\,\partial_y N_j)\). Zeilensummen null, wie bei der Feder.
Assemblieren, Randwerte, lösen liefert für den symmetrischen Träger die Mittelspalte auf 160 °C: 2D ist hier noch 1D in Verkleidung.

Roter Faden

Zurück: Der ganze Rechenweg ist die Wiederkehr von Kapitel 7 — nur der Interpolationsbaustein wechselte von Streckenanteilen (Hütchen) zu Flächenanteilen (Zelt). Die Elementmatrix aus Steigung mal Steigung, die Assemblierung durch Aufaddieren, das Streich-Rezept für Randwerte (Kapitel 8) und der Gauß-Seidel-Löser (Kapitel 4) sind alle unverändert übernommen. Selbst die Probe „Zeilensumme null” ist dieselbe wie bei der Feder (Kapitel 3).

Vor: In Kapitel 11 lösen Rechtecke die Dreiecke ab — mit bilinearer statt linearer Interpolation, und wir vergleichen beide. Kapitel 12 fragt, wie genau diese groben Netze eigentlich sind (Netzkonvergenz), denn hier war die Lösung nur darum exakt, weil sie zufällig eine Gerade war. Und dasselbe Netz, das hier die Temperatur trug, trägt in Kapitel 14 die Verschiebungen des sich verbiegenden Trägers — dann sitzen an jedem Knoten zwei Zahlen statt einer.

Übungen

Ü 10.1 (Verstehen). Gegeben das Dreieck mit den Ecken \(V_1 = (0,0)\), \(V_2 = (6,0)\), \(V_3 = (0,3)\) und der Punkt \(P = (2,1)\). Berechne von Hand die drei Flächenanteile (halbe Kreuzdifferenz) und prüfe die Summe. Kontrolliere danach im Dreiecks-Labor (Teil 1).

Musterlösung zu Ü 10.1

Gesamtfläche: \(A = \tfrac12\,|(6-0)(3-0) - (0-0)(0-0)| = \tfrac12\cdot 18 = 9\).

\(A_1 = \text{Fläche}(P, V_2, V_3) = \tfrac12\,|(6-2)(3-1) - (0-2)(0-1)| = \tfrac12\,|8 - 2| = 3\) → \(\lambda_1 = 3/9 = \tfrac13\).
\(A_2 = \text{Fläche}(P, V_3, V_1) = \tfrac12\,|(0-2)(0-1) - (0-2)(3-1)| = \tfrac12\,|2 + 4| = 3\) → \(\lambda_2 = 3/9 = \tfrac13\).
\(A_3 = \text{Fläche}(P, V_1, V_2) = \tfrac12\,|(0-2)(0-1) - (6-2)(0-1)| = \tfrac12\,|2 + 4| = 3\) → \(\lambda_3 = 3/9 = \tfrac13\).

Probe: \(\tfrac13 + \tfrac13 + \tfrac13 = 1\). ✓ Der Punkt \(P = (2,1)\) ist der Schwerpunkt des Dreiecks (Mittel der drei Ecken: \(((0+6+0)/3,\ (0+0+3)/3) = (2,1)\)) — deshalb alle Anteile gleich.

Ü 10.2 (Verändern). Jetzt kommt die Wärme nicht mehr von der Seite, sondern von unten — die Brandszene des Buches. Beflamme die Unterkante des Trägers: Halte die untere Kante (Knoten 1, 2, 3) auf 300 °C (die Brandgase) und die obere Kante (Knoten 7, 8, 9) auf 20 °C (die kühle Betondecke drückt von oben); die Mittelreihe (Knoten 4, 5, 6) bleibt frei. Sage vorher: Bleibt die Temperatur über die Höhe gleich wie im Kapitel, oder entsteht ein Gefälle von unten nach oben? Löse und deute das Farbbild. Das Skript loesungen/kap10_ue2.py rechnet es vor.

Musterlösung zu Ü 10.2

Jetzt sind Unter- und Oberkante festgehalten, die Mittelreihe ist frei. Aus der Symmetrie stellt sich die Mittelreihe genau dazwischen ein: alle drei Knoten (4, 5, 6) auf 160 °C. Das Feld ist dasselbe wie im Kapitel — nur um eine Vierteldrehung gekippt: Vorher war die Temperatur auf jeder senkrechten Linie konstant und wuchs über die Länge; jetzt ist sie auf jeder waagerechten Linie konstant und fällt über die Höhe, von 300 °C an der beflammten Unterkante auf 20 °C an der kühlen Oberkante.

Und genau das ist die Vorschau auf das Finale. Ein solches Höhengefälle — heiße Unterseite, kühle Oberseite — dehnt den unteren Rand stärker als den oberen und krümmt den Träger nach oben (thermisches Verbiegen, Kapitel 14, Kapitel 15). Die Wärmerechnung liefert hier das Gefälle; die Verformungsrechnung macht später die Biegung daraus. Das Feld, auf dem der Träger sich verbiegt, steht damit schon.

Ü 10.3 (Übertragen). Schreibe ein eigenes L-förmiges Mini-Netz aus vier Dreiecken als Knoten- und Dreiecksliste auf (ein \(2\times 2\)-Quadrat ohne die rechte obere Zelle — sechs Ecken, vier Dreiecke) und schicke es durch die Maschinerie: linke Kante 20 °C, rechte Kante 300 °C. Wo stellt sich die einspringende Ecke ein? Das Skript loesungen/kap10_ue3.py zeigt eine Musterlösung.

Musterlösung zu Ü 10.3

Ein möglicher Bauplan (gegen den Uhrzeigersinn nummeriert):

Knoten	1	2	3	4	5	6
\((x, y)\)	(0,0)	(2,0)	(2,1)	(1,1)	(1,2)	(0,2)

Dreiecke: \(1\text{–}2\text{–}4\), \(1\text{–}4\text{–}6\), \(2\text{–}3\text{–}4\), \(4\text{–}5\text{–}6\) (Flächen \(1, 1, \tfrac12, \tfrac12\), Summe 3 — das L aus drei Einheitsquadraten). Hält man die linke Kante (Knoten 1, 6) auf 20 °C und die rechte (Knoten 2, 3) auf 300 °C, stellt sich die einspringende Ecke (Knoten 4) auf rund 122 °C ein — näher an der kalten Seite, weil sie an das große kalte Gebiet links grenzt. Die entscheidende Lehre: Der einzige Aufwand war, die zwei Listen richtig aufzuschreiben; die Maschinerie lief unverändert.

Das Kleingedruckte

Drei ehrliche Feinheiten.

Erstens — dünne Dreiecke. Unsere Träger-Dreiecke sind mit \(0{,}5\ \mathrm{m} \times 0{,}05\ \mathrm{m}\) ziemlich langgestreckt (Verhältnis 10 zu 1). Das sieht man den großen Zahlen (25,25, 25) in der Elementmatrix an: Die steilen \(y\)-Steigungen entstehen aus der kurzen Höhe. Für unser gutmütiges, gleichmäßiges Problem ist das kein Problem, aber als Faustregel gilt: Vermeide Nadeln. Sehr spitze Dreiecke verschlechtern die Genauigkeit und machen das Gleichungssystem schwerer lösbar. Netzgeneratoren (Delaunay) achten deshalb auf gedrungene Dreiecke. Wie man Netzqualität misst, gehört zu Kapitel 12.

Zweitens — nur lineare Elemente. Wir benutzen ausschließlich P1-Dreiecke mit gerader Interpolation. Über einem Dreieck ist das Feld deshalb immer eine Ebene, sein Gradient konstant. Elemente höherer Ordnung (gekrümmte Ansätze) können mit weniger Dreiecken mehr, sind aber ein Thema für sich und kommen in diesem Buch nur als Ausblick vor.

Drittens — die Exaktheit ist wieder geliehen. Dass die Mittelspalte exakt 160 °C traf, liegt wie in Kapitel 7 daran, dass die wahre Lösung selbst eine Ebene ist (linear in \(x\), unabhängig von \(y\)) — und eine Ebene stellen die P1-Dreiecke exakt dar. Sobald die Wahrheit krumm ist (eine einspringende Ecke, ein Loch, ein runder Rand), trifft die FEM nur noch näherungsweise, und man braucht mehr Dreiecke. Wie viel genauer ein feineres Netz wird — und wie man das ehrlich misst —, ist das Thema von Kapitel 12.