9  Jetzt mit Zeit

Bisher stand unser Träger immer schon fertig aufgeheizt vor uns. Wir haben die Gerade berechnet, auf der sich die Temperaturen im eingeschwungenen Zustand einpendeln — 20 °C am kühlen Auflager links, 300 °C beim Brandherd rechts, sauber dazwischen gestaffelt. Aber die entscheidende Frage des Brandschutzes haben wir dabei übersprungen: Wie lange dauert es überhaupt, bis es so weit ist? Wie viele Minuten kriecht die Hitze nach links, ehe sie auch das ferne Ende im Auflager erreicht?

Ab jetzt tickt die Uhr. In Kapitel 1 lief diese Simulation schon einmal vor deinen Augen ab — die Wärmefront, die vom Brandherd ins Stahlprofil kriecht —, aber ihr Maschinenraum war eine Blackbox. Dieses Kapitel baut genau diesen Maschinenraum. Und es hält ein Versprechen ein, das dort schon im Code steckte: die kleine Schleife mit den „expliziten Zeitschritten”. Am Ende dieses Kapitels weißt du, warum sie so aussieht, wie sie aussieht — und warum sie, einen Tick zu mutig eingestellt, mit einem Knall explodiert.

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

lam, A, L = 50.0, 0.01, 1.0
rho, c = 7850.0, 460.0
knoten = 21
h = L / (knoten - 1)
k = lam * A / h
C = rho * c * A * h            # innere Knotenkapazitaet (gelumpt)
dt = 600.0

# Implizite Schritte (immer stabil) -- reine Bandmatrix pro Schritt.
haupt = (C / dt + 2 * k) * np.ones(knoten - 2)
neben = -k * np.ones(knoten - 3)
M = np.diag(haupt) + np.diag(neben, 1) + np.diag(neben, -1)

T = np.full(knoten, 20.0)
T[-1] = 300.0
schritte = 60
verlauf, zeiten = [T.copy()], [0.0]
for s in range(schritte):
    rs = (C / dt) * T[1:-1].copy()
    rs[0] += k * 20.0
    rs[-1] += k * 300.0
    T[1:-1] = np.linalg.solve(M, rs)
    verlauf.append(T.copy())
    zeiten.append((s + 1) * dt)

orte = [i * h for i in range(knoten)]
gerade = [20.0 + (i / (knoten - 1)) * 280.0 for i in range(knoten)]

fig, ax = plt.subplots(figsize=(8.0, 4.0))

def zeichne(frame):
    ax.clear()
    ax.plot(orte, gerade, "--", color="gray", lw=1.3,
            label="eingeschwungene Gerade")
    ax.plot(orte, verlauf[frame], "o-", color="tab:red", ms=6,
            markeredgecolor="black", lw=2.0)
    ax.text(0.05, 305, "nach %4.1f Stunden" % (zeiten[frame] / 3600.0),
            fontsize=12, va="top")
    ax.set_xlim(-0.03, 1.03)
    ax.set_ylim(0, 330)
    ax.set_xlabel("Ort entlang des Trägers (m)")
    ax.set_ylabel("Temperatur (°C)")
    ax.set_title("Die Wärmefront wandert nach links")
    ax.legend(loc="center right")
    return []

ani = FuncAnimation(fig, zeichne, frames=len(verlauf), interval=140,
                    blit=False)
plt.close(fig)
HTML(ani.to_jshtml())
Abbildung 9.1: Die Wärmefront wandert: Start ist der kalte Träger mit nur dem rechten Ende beim Brandherd (300 °C). Über rund zehn Stunden schiebt sich die Wärme nach links, bis die eingeschwungene Gerade (gestrichelt) erreicht ist. Zeitraffer: zehn Stunden in wenigen Sekunden; die eingeblendete Uhr zählt die echten Stunden.

Was man hier sieht: einen Träger, der nicht auf einen Schlag heiß ist, sondern Minute um Minute wärmer wird. Anfangs steht nur das rechte Ende beim Brandherd, der Rest liegt kalt im Auflager. Dann frisst sich die Wärme nach links, erst steil, dann immer träger, bis sich das Profil an die vertraute Gerade aus Kapitel 6 anschmiegt und stehen bleibt. Wichtig: Nach 30 und 90 Minuten — den ersten beiden Bildern der Reihe — ist erst das rechte Drittel warm; das Auflager bleibt lange kalt. Diese ganze Bewegung — den Weg zur Geraden — wollen wir jetzt selbst rechnen.

Lernziele

Nach diesem Kapitel kannst du …

  1. … den Speicherterm \(\rho c\,\partial T/\partial t\) physikalisch erklären (die Bilanz aus Kapitel 6, bei der das Konto sich jetzt doch ändert) und die Rolle der volumetrischen Wärmekapazität \(\rho c\) benennen,
  2. … erklären, was die Massenmatrix speichert und was Lumping damit tut — die Wärmekapazität eines Elements auf seine beiden Knoten verteilen,
  3. … einen expliziten Euler-Schritt von Hand rechnen — neuer Wert = alter Wert + Zeitschritt · momentane Änderungsrate — am 5-Knoten-Träger mit echten Zahlen,
  4. … das Stabilitätslimit erklären und als Faustregel formulieren: halbe Maschenweite erzwingt ein Viertel des Zeitschritts, sonst explodiert die Rechnung,
  5. … den impliziten Schritt erklären — alle neuen Werte gleichzeitig aus einem Gleichungssystem bestimmen — und begründen, wann sich sein Mehraufwand lohnt.
WarnungNaheliegende Vermutung

Vermutung: „Der Zeitschritt ist wie die Netzfeinheit nur eine Frage der Genauigkeit. Wähle ich ihn grob, wird das Ergebnis eben etwas ungenauer — aber es bleibt ein Ergebnis.”

Warum sie naheliegt: Bei allem, was wir bisher genähert haben — die Fläche unter einer Kurve, das Temperaturprofil zwischen den Knoten —, galt genau das: grob = ungenau, fein = genau. Ein sanfter Regler, kein Abgrund.

Was stattdessen stimmt: Beim expliziten Zeitschritt gibt es eine scharfe Grenze. Diesseits davon wird ein größerer Schritt tatsächlich nur etwas ungenauer. Aber einen Tick jenseits der Grenze wird die Rechnung nicht ungenau, sondern sie explodiert: Die Temperaturen fangen an zu zappeln, das Zappeln schaukelt sich von Schritt zu Schritt auf, und nach wenigen Sekunden stehen an den Knoten Millionen Grad — mal plus, mal minus. Ein zu großer Schritt überzieht ein Konto, und die Überziehung wächst sich zur Katastrophe aus. Das rechnen und erleben wir in diesem Kapitel.

9.1 Das Konto ändert sich jetzt doch

In Kapitel 6 war ein Trägerstück ein Girokonto für Wärme: vom Brandherd strömen 140 W herein, zum kühlen Auflager fließen 140 W hinaus, die Bilanz geht auf null auf, und die Temperatur des Stücks bleibt stehen. Das war die Bedingung des eingeschwungenen Zustands — und genau deshalb galt sie erst am Ende des Aufheizens.

Solange der Träger noch kalt ist, stimmt die Bilanz eben nicht. An einem Stück nahe dem Brandherd strömt mehr Wärme herein, als nach links weiterfließt. Wohin geht der Überschuss? Er bleibt im Stück und heizt es auf. Das ist keine Nebensache, sondern der ganze Vorgang: Die Differenz zwischen Zufluss und Abfluss ist genau das, was die Temperatur des Stücks von Sekunde zu Sekunde ändert.

Code
import matplotlib.pyplot as plt

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

# Das freigeschnittene Stueck (leicht ins Warme verfaerbt).
ax.add_patch(plt.Rectangle((0.25, 0.0), 0.5, 1.0, facecolor="#fbcda6",
                           edgecolor="black", lw=1.5))
ax.text(0.5, 0.42, "Trägerstück\num den Knoten\nbei 0,5 m", ha="center",
        va="center", fontsize=11)

for xk in [0.25, 0.5, 0.75]:
    ax.plot(xk, 0.0, "o", color="tab:blue", ms=10, markeredgecolor="black",
            clip_on=False, zorder=6)

# Waermestrom rein (vom Brandherd rechts): 140 W
ax.annotate("", xy=(0.80, 0.72), xytext=(1.05, 0.72),
            arrowprops=dict(arrowstyle="->", color="tab:red", lw=3.0))
ax.text(0.92, 0.86, "rein: 140 W", color="tab:red", ha="center", fontsize=12)

# Waermestrom raus (nach links zum Auflager): nur 100 W
ax.annotate("", xy=(-0.05, 0.72), xytext=(0.20, 0.72),
            arrowprops=dict(arrowstyle="->", color="tab:red", lw=2.4))
ax.text(0.075, 0.86, "raus: 100 W", color="tab:red", ha="center", fontsize=12)

# Der neue Speicherpfeil (nach oben): die Differenz laedt das Stueck auf.
ax.annotate("", xy=(0.5, 1.28), xytext=(0.5, 0.55),
            arrowprops=dict(arrowstyle="-|>", color="tab:purple", lw=3.0))
ax.text(0.62, 1.05, "Speicher:\n40 W heizen auf", color="tab:purple",
        ha="left", va="center", fontsize=11)

ax.text(0.5, -0.28, "Bilanz:  rein − raus  =  140 − 100  =  40 W  ≠  0",
        ha="center", fontsize=13, color="black",
        bbox=dict(boxstyle="round", facecolor="#f0e6f5",
                  edgecolor="tab:purple"))

ax.set_xlim(-0.15, 1.30)
ax.set_ylim(-0.45, 1.45)
ax.axis("off")
plt.tight_layout()
plt.show()
Abbildung 9.2: Dieselbe Bilanz wie in Kapitel 6, aber im Aufheizen: Vom Brandherd strömen 140 W herein, nach links zum kühlen Auflager fließen erst 100 W ab. Die Differenz von 40 W verschwindet nicht — sie lädt das Stück auf (der lila Pfeil nach oben) und hebt seine Temperatur. Dieser Speicherpfeil ist der neue Term dieses Kapitels.

Was man hier sieht: dasselbe Konto wie in Kapitel 6, aber mitten im Aufheizen. Diesmal ist der Abfluss kleiner als der Zufluss, und die Differenz — hier 40 W — ist nicht verloren. Sie steckt als gespeicherte Wärme im Stück und treibt seine Temperatur nach oben (der lila Pfeil). Genau diesen Pfeil hat Kapitel 6 bewusst weggelassen, um „ein schweres Ding nach dem anderen” zu behandeln. Jetzt holen wir ihn.

9.1.1 Wie viel Wärme steckt in einem Grad?

Wie stark ein bisschen überschüssige Wärme die Temperatur hebt, hängt vom Material ab. Zwei Stoffgrößen des Baustahls bestimmen das: die Dichte \(\rho = 7850\ \mathrm{kg/m^3}\) (wie viel Masse in einem Kubikmeter steckt) und die spezifische Wärmekapazität \(c = 460\ \mathrm{J/(kg\,K)}\) (wie viel Wärme ein Kilogramm um ein Grad erwärmt). Multipliziert man beide, bekommt man die volumetrische Wärmekapazität:

\[ \rho\,c \;=\; 7850 \cdot 460 \;=\; 3\,611\,000\ \frac{\mathrm{J}}{\mathrm{m^3\,K}} \;=\; 3{,}61\ \frac{\mathrm{MJ}}{\mathrm{m^3\,K}} . \]

Lies das so: Um einen Kubikmeter Stahl um ein einziges Grad zu erwärmen, muss man 3,61 Millionen Joule hineinstecken. Das ist der Preis der Trägheit — und der Grund, warum der Träger sich Zeit lässt, während der Brand unter ihm längst lodert.

In der Sprache der Wärmeleitungsgleichung tritt dieser Speicher als neuer Term links neben die alte Bilanz. Statt „Zufluss minus Abfluss = 0” heißt es jetzt „Zufluss minus Abfluss = das, was gespeichert wird”:

\[ \underbrace{\rho\,c\,\frac{\partial T}{\partial t}}_{\text{Speicher (neu)}} \;=\; \underbrace{\frac{\partial}{\partial x} \left( \lambda\,\frac{\partial T}{\partial x} \right)}_{\text{Zufluss − Abfluss (Kapitel 6)}} . \]

Das Zeichen \(\partial T/\partial t\) ist die Änderungsrate der Temperatur über der Zeit — Grad pro Sekunde, genau die Steigung aus Kapitel 5, nur dass die Waagerechte jetzt nicht der Ort \(x\) ist, sondern die Zeit \(t\). Rechts steht buchstäblich die Bilanz aus Kapitel 6. Im eingeschwungenen Zustand ist die linke Seite null (nichts ändert sich mehr), und wir sind zurück bei der alten Geraden. Alles Neue an diesem Kapitel steckt in dem einen linken Term.

9.2 Von der Wärmekapazität zur Massenmatrix

Den Zufluss-minus-Abfluss-Term haben wir in Kapitel 7 in die Maschinerie übersetzt: Er wurde zur Steifigkeitsmatrix mit der Wärme-Steifigkeit \(k = \lambda A/h = 2\ \mathrm{W/K}\) zwischen je zwei Knoten. Mit dem neuen Speicherterm machen wir jetzt genau dasselbe Spiel — dieselbe schwache Form, dieselben Hütchenfunktionen \(N_i\) aus Kapitel 7.

Setzt man den Ansatz \(T(x,t) = \sum_j T_j(t)\,N_j(x)\) (Knotenwerte mal Hütchen, jetzt mit der Zeit veränderliche Knotenwerte) in den Speicherterm ein und gewichtet wie immer mit dem Hütchen \(N_i\), entsteht aus \(\int \rho c\, N_i\, \partial T/\partial t\, \mathrm{d}x\) eine Summe: Jeder Knoten \(j\) trägt mit seiner Änderungsrate \(\mathrm{d}T_j/\mathrm{d}t\) bei, gewichtet mit dem Überlapp der beiden Hütchen. Diese Gewichte bilden eine neue Matrix — die Massenmatrix \(\mathbf{M}\). Für ein einzelnes Element sieht sie so aus:

TippErgebnis-Kasten: die konsistente Elementmassenmatrix

\[ \mathbf{M}_{\text{Element}} \;=\; \frac{\rho\,c\,A\,h}{6} \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} . \]

Sie wiegt ab, wie stark die Temperaturänderung an einem Knoten in die Bilanz seines Nachbarn hineinwirkt. Weil sich benachbarte Hütchen überlappen, sind auch die Nebeneinträge (die \(1\)en) ungleich null: Die Änderungsraten benachbarter Knoten sind gekoppelt. Diese exakt integrierte Form heißt konsistente Massenmatrix. Wo die \(6\) und die \(2, 1\) herkommen, steht im Kleingedruckten am Kapitelende — für den Weg durch dieses Kapitel brauchen wir sie nicht.

Der Vorfaktor ist die gesamte Wärmekapazität eines Elements: Masse mal spezifische Wärme, also \(\rho\,c\,A\,h\). Mit den Zahlen unseres 5-Knoten-Trägers (\(A = 0{,}01\ \mathrm{m^2}\), \(h = 0{,}25\ \mathrm{m}\)):

\[ \rho\,c\,A\,h \;=\; 3\,611\,000 \cdot 0{,}01 \cdot 0{,}25 \;=\; 9027{,}5\ \frac{\mathrm{J}}{\mathrm{K}} . \]

Ein Elementstück des Trägers braucht also 9027,5 Joule, um ein Grad wärmer zu werden. Diese Zahl trägt das ganze Kapitel.

9.2.1 Lumping: jeder Knoten bekommt seinen Eimer

Die konsistente Massenmatrix hat einen Haken: Ihre Nebeneinträge koppeln die Knoten. Wollen wir gleich einen Zeitschritt machen, müssten wir wegen dieser Kopplung selbst bei der einfachsten Rechnung ein Gleichungssystem lösen — der Vorteil eines „billigen” expliziten Schritts wäre dahin.

Es gibt einen ehrlichen Ausweg, und er hat einen anschaulichen Namen: Lumping, vom englischen lump, der Klumpen. Statt die Speicherwirkung über das Element zu verschmieren, klumpen wir sie an den Knoten zusammen. Jedes Element gibt seine Kapazität \(\rho c A h\) je zur Hälfte an seine beiden Knoten ab — genau die Halbe-Halbe-Logik, mit der schon in Kapitel 7 die Wärmequelle auf die Knoten verteilt wurde. Das Bild dazu: An jedem Knoten hängt ein Eimer, in dem sich seine Wärme sammelt.

Code
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10.0, 3.8))

orte = [0.0, 0.25, 0.5, 0.75, 1.0]
kapazitaeten = [4513.75, 9027.5, 9027.5, 9027.5, 4513.75]

# Der Traeger als Linie mit Knoten und Elementgrenzen.
ax.plot([0, 1], [0, 0], "-", color="black", lw=2.0, zorder=1)
for x in orte:
    ax.plot(x, 0, "o", color="tab:blue", ms=12, markeredgecolor="black",
            zorder=3)

# Ein „Eimer" (Rechteck) je Knoten, Fuellhoehe proportional zur Kapazitaet.
for x, kap in zip(orte, kapazitaeten):
    hoehe = kap / 9027.5
    breite = 0.10
    ax.add_patch(plt.Rectangle((x - breite / 2, 0.25), breite, hoehe,
                               facecolor="#9ecae1", edgecolor="black",
                               lw=1.4, zorder=2))
    ax.text(x, 0.25 + hoehe + 0.06, "%.1f\nJ/K" % kap, ha="center",
            va="bottom", fontsize=9)

# Elementbeschriftung.
for e in range(4):
    xm = (orte[e] + orte[e + 1]) / 2
    ax.text(xm, -0.18, "Element %d\n9027,5 J/K" % (e + 1), ha="center",
            va="top", fontsize=8, color="gray")

ax.set_xlim(-0.1, 1.1)
ax.set_ylim(-0.45, 1.55)
ax.axis("off")
ax.set_title("Die Wärmekapazität, auf Knoten-Eimer verteilt")
plt.tight_layout()
plt.show()
Abbildung 9.3: Lumping am 5-Knoten-Träger: Jedes der vier Elemente (Kapazität 9027,5 J/K) gibt die Hälfte an jeden seiner beiden Knoten ab. Innere Knoten bekommen von zwei Elementen je 4513,75 J/K — zusammen ein voller Eimer von 9027,5 J/K. Die beiden Randknoten grenzen nur an ein Element und tragen den halben Eimer, 4513,75 J/K.

Was man hier sieht: fünf Eimer, aber nicht gleich groß. Die drei inneren Knoten grenzen an zwei Elemente und sammeln von jedem die Hälfte ein — \(4513{,}75 + 4513{,}75 = 9027{,}5\ \mathrm{J/K}\), ein voller Eimer. Die beiden Randknoten sitzen am Ende und grenzen nur an ein einziges Element; ihr Eimer ist halb so groß, \(4513{,}75\ \mathrm{J/K}\). Wichtig: Die Gesamtsumme bleibt erhalten — nichts geht verloren, die Wärme wird nur anders verteilt. Aus der Massenmatrix ist eine reine Diagonale geworden, ein Eimer je Knoten, keine Kopplung mehr. Genau das brauchen wir für den nächsten Schritt.

9.3 Der explizite Euler-Schritt

Jetzt haben wir alle Teile beisammen. Die Bilanz an einem inneren Knoten lautet in Worten: Die gespeicherte Wärme ändert sich so schnell, wie netto Wärme hineinfließt. In Zahlen — Eimer mal Änderungsrate = Zufluss minus Abfluss:

\[ C_i \cdot \frac{\mathrm{d}T_i}{\mathrm{d}t} \;=\; k\,\bigl(T_{i-1} - 2\,T_i + T_{i+1}\bigr) . \]

Links steht der Eimer \(C_i\) (für innere Knoten \(9027{,}5\ \mathrm{J/K}\)) mal der Änderungsrate. Rechts steht die Bilanz aus Kapitel 6 und Kapitel 7, ausgedrückt durch die drei Knoten: Der Knoten schaut nach links und rechts, und was seine beiden Nachbarn zusammen mehr haben als das Doppelte von ihm, treibt ihn. Das Minus-Doppelte ist genau die \(-2, +4, -2\)-Zeile der Steifigkeitsmatrix.

Der Trick, aus dieser Beziehung eine Vorhersage zu machen, ist der explizite Euler-Schritt, benannt nach Leonhard Euler (Euler 1768). Er ist bestechend einfach: Ersetze die Änderungsrate durch „Wie viel ändert sich pro Zeitschritt \(\Delta t\)” und löse nach dem neuen Wert auf:

\[ T_i^{\text{neu}} \;=\; T_i^{\text{alt}} \;+\; \frac{\Delta t}{C_i}\; k\,\bigl(T_{i-1} - 2\,T_i + T_{i+1}\bigr) . \]

In Worten, und das ist der ganze Satz, den man sich merken muss:

Neuer Wert = alter Wert + Zeitschritt · momentane Änderungsrate.

Jeder Knoten schaut auf seine zwei alten Nachbarwerte, rechnet seine Änderungsrate aus und macht einen Schritt in diese Richtung. Kein Gleichungssystem, keine Kopplung — nur eine Division je Knoten. Genau das macht das Lumping so wertvoll.

WichtigVorhersage-Punkt

Bevor du weiterliest: Gleich rechnen wir den Träger von Hand los. Der Startzustand ist die Liste [20, 20, 20, 20, 300] aus Kapitel 2 — alles kalt, nur das rechte Ende beim Brandherd. Wir wählen \(\Delta t = 600\ \mathrm{s}\) (zehn Minuten). Welcher der drei inneren Knoten (2, 3 oder 4) ändert sich im allerersten Schritt als Einziger? Und warum die anderen beiden (noch) nicht? Leg dich fest.

9.3.1 Drei Schritte von Hand

Wir rechnen mit \(\Delta t = 600\ \mathrm{s}\), \(k = 2\ \mathrm{W/K}\) und dem inneren Eimer \(C = 9027{,}5\ \mathrm{J/K}\). Der Vorfaktor ist für alle inneren Knoten gleich:

\[ \frac{\Delta t}{C}\,k \;=\; \frac{600}{9027{,}5}\cdot 2 \;\approx\; 0{,}1329 . \]

Die Randknoten 1 (20 °C, Auflager) und 5 (300 °C, Brandherd) sind festgehalten — Dirichlet, wie in Kapitel 8 — und ändern sich nie. Nur die Knoten 2, 3 und 4 rechnen. Nehmen wir den ersten Schritt für Knoten 4 (rechts, neben dem Brandherd):

\[ T_4^{\text{neu}} = 20 + 0{,}1329\cdot\bigl(\underbrace{20}_{T_3} - 2\cdot 20 + \underbrace{300}_{T_5}\bigr) = 20 + 0{,}1329\cdot 280 = 57{,}2\ °\mathrm{C} . \]

Für Knoten 2 und 3 dagegen sind im ersten Schritt beide Nachbarn noch bei 20 °C, die Klammer wird \(20 - 40 + 20 = 0\) — sie ändern sich nicht. Die Wärme muss erst bei ihnen ankommen. Genau das war die Vorhersage-Frage.

Code
import matplotlib.pyplot as plt

zeilen = [
    ("Start (0 s)",   [20.0, 20.0, 20.0, 20.0, 300.0]),
    ("nach 600 s",    [20.0, 20.0, 20.0, 57.220, 300.0]),
    ("nach 1200 s",   [20.0, 20.0, 24.947, 84.544, 300.0]),
    ("nach 1800 s",   [20.0, 20.658, 32.212, 105.262, 300.0]),
]
# Welche inneren Knoten haben sich gegenueber der Vorzeile veraendert?
neu_bewegt = [
    [False, False, False, False, False],
    [False, False, False, True,  False],
    [False, False, True,  True,  False],
    [False, True,  True,  True,  False],
]

fig, ax = plt.subplots(figsize=(9.5, 3.2))
ax.axis("off")
spalten = ["Zeitpunkt", "Knoten 1", "Knoten 2", "Knoten 3", "Knoten 4",
           "Knoten 5"]
tabelle = ax.table(
    cellText=[[name] + ["%.1f" % w for w in werte]
              for name, werte in zeilen],
    colLabels=spalten, loc="center", cellLoc="center")
tabelle.auto_set_font_size(False)
tabelle.set_fontsize(11)
tabelle.scale(1.0, 1.7)

for r, werte in enumerate(zeilen):
    for c in range(5):
        zelle = tabelle[r + 1, c + 1]
        if neu_bewegt[r][c]:
            zelle.set_facecolor("#ffd6d6")
        if c == 0 or c == 4:
            zelle.set_facecolor("#e0e0e0")   # feste Randknoten grau
ax.set_title("Explizite Handrechnung: die Wärmefront kriecht nach links")
plt.tight_layout()
plt.show()
Abbildung 9.4: Die ersten drei expliziten Zeitschritte von Hand, Δt = 600 s. Rot markiert ist in jedem Schritt, welche Knoten sich neu bewegt haben. Die Wärme kriecht vom Knoten 4 (am Brandherd) Schritt für Schritt nach links — nach einer halben Stunde hat sie gerade eben Knoten 2 erreicht.

Was man hier sieht: die Front als Zahlenwelle. Im ersten Schritt bewegt sich nur Knoten 4. Im zweiten hat die Wärme Knoten 3 erreicht, im dritten auch Knoten 2. Die beiden grauen Randspalten stehen fest. Jede Zahl in dieser Tabelle ist mit dem Taschenrechner nachprüfbar — es steckt kein einziger versteckter Schritt darin.

Und genau diese Tabelle ist eine for-Schleife. Der Kern des expliziten Schritts sind acht Zeilen reines Python — eine je Nachbar-Blick, so wie in der Handrechnung:

Code
def expliziter_schritt(temperaturen, kapazitaeten, k, dt):
    """Ein expliziter Euler-Zeitschritt am 1D-Träger.

    temperaturen: Liste der jetzigen Knotentemperaturen (Grad Celsius).
    kapazitaeten: gelumpte Knotenkapazitaeten C_i (J/K).
    k: Waerme-Steifigkeit lambda*A/h (W/K).
    dt: Zeitschritt (Sekunden).
    Rueckgabe: neue Temperaturliste; die Randknoten bleiben fest.
    """
    neu = list(temperaturen)
    for i in range(1, len(temperaturen) - 1):
        zufluss = k * (temperaturen[i - 1]
                       - 2.0 * temperaturen[i]
                       + temperaturen[i + 1])
        neu[i] = temperaturen[i] + dt / kapazitaeten[i] * zufluss
    return neu


# Die Handrechnung nachfahren: Start [20, 20, 20, 20, 300], drei Schritte.
kapazitaeten = [4513.75, 9027.5, 9027.5, 9027.5, 4513.75]
temperaturen = [20.0, 20.0, 20.0, 20.0, 300.0]
for schritt in range(3):
    temperaturen = expliziter_schritt(temperaturen, kapazitaeten, 2.0, 600.0)
    print("nach %4d s:" % ((schritt + 1) * 600), end=" ")
    for wert in temperaturen:
        print("%8.3f" % wert, end="")
    print()
nach  600 s:   20.000  20.000  20.000  57.220 300.000
nach 1200 s:   20.000  20.000  24.947  84.544 300.000
nach 1800 s:   20.000  20.658  32.212 105.262 300.000

Interpretation: Zeile für Zeile dieselben Zahlen wie in der Tabelle — 57.220, dann 24.947 und 84.544, dann 20.658, 32.212 und 105.262. Die for-Schleife über i ist der Blick jedes Knotens auf seine Nachbarn i-1 und i+1; die Division durch kapazitaeten[i] ist der Eimer. Mehr passiert bei einem expliziten Zeitschritt nicht.

9.4 Das Experiment: den Zeitschritt vergrößern

Zehn Minuten pro Handrechnung sind mühsam. Bis der Träger eingeschwungen ist, vergehen viele Stunden (dazu gleich mehr) — das sind Tausende von 600-Sekunden-Schritten. Der naheliegende Wunsch: den Zeitschritt größer machen, dann sind es weniger Schritte. Probieren wir es. Was passiert bei \(\Delta t = 1800\ \mathrm{s}\) (einer halben Stunde), und was bei \(\Delta t = 3000\ \mathrm{s}\) (fünfzig Minuten)?

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

k = 2.0
kap = np.array([4513.75, 9027.5, 9027.5, 9027.5, 4513.75])
orte = [0.0, 0.25, 0.5, 0.75, 1.0]

def lauf(dt, schritte):
    T = np.array([20.0, 20.0, 20.0, 20.0, 300.0])
    verlauf = [T.copy()]
    for s in range(schritte):
        neu = T.copy()
        for i in [1, 2, 3]:
            neu[i] = T[i] + dt / kap[i] * k * (T[i-1] - 2*T[i] + T[i+1])
        T = neu
        verlauf.append(T.copy())
    return verlauf

stabil = lauf(1800.0, 18)
explo = lauf(3000.0, 18)

fig, (axl, axr) = plt.subplots(1, 2, figsize=(11.0, 4.2))

def zeichne(frame):
    axl.clear(); axr.clear()
    axl.axhspan(20, 300, color="#e8f5e9", zorder=0)
    axl.plot(orte, stabil[frame], "o-", color="tab:green", ms=8,
             markeredgecolor="black")
    axl.set_ylim(0, 320)
    axl.set_title("Δt = 1800 s  (stabil)")
    axl.set_xlabel("Ort (m)"); axl.set_ylabel("Temperatur (°C)")

    axr.axhspan(20, 300, color="#fdecea", zorder=0)
    axr.plot(orte, explo[frame], "o-", color="tab:red", ms=8,
             markeredgecolor="black")
    grenze = max(400.0, np.max(np.abs(explo[frame])) * 1.15)
    axr.set_ylim(-grenze, grenze)
    axr.set_title("Δt = 3000 s  (explodiert)")
    axr.set_xlabel("Ort (m)"); axr.set_ylabel("Temperatur (°C)")
    fig.suptitle("Schritt %d" % frame, fontsize=13)
    return []

ani = FuncAnimation(fig, zeichne, frames=len(stabil), interval=380,
                    blit=False)
plt.close(fig)
HTML(ani.to_jshtml())
Abbildung 9.5: Zwei Läufe desselben 5-Knoten-Trägers, nebeneinander. Links Δt = 1800 s: die Front wandert brav ein. Rechts Δt = 3000 s: erst zappeln die Werte nur, dann schießen sie über die 300 °C hinaus, kippen ins Negative und sprengen nach wenigen Schritten die Skala (beachte die y-Achse). Das ist die Explosion.

Was man hier sieht: zwei Rechnungen mit demselben Programm, demselben Träger, demselben Startzustand — nur der Zeitschritt ist anders. Links mit \(\Delta t = 1800\ \mathrm{s}\) läuft alles ruhig ein. Rechts mit \(\Delta t = 3000\ \mathrm{s}\) beginnt es harmlos, dann fangen die inneren Knoten an, gegeneinander zu schwingen: einer schießt über 300 °C hinaus, sein Nachbar fällt unter null, im nächsten Schritt ist es umgekehrt, nur heftiger. Nach wenigen Schritten stehen dort Hunderte, bald Tausende Grad — mit wechselndem Vorzeichen, und es wächst weiter. Physikalisch ist das Unsinn: Nichts am Träger kann kälter als das Auflager oder heißer als der Brandherd werden. Es ist kein Messfehler und keine Ungenauigkeit, sondern ein Versagen des Verfahrens.

HinweisZappeln ist nicht gleich Zappeln

In Kapitel 4 hat die Lösung beim Gauß-Seidel-Verfahren auch „gezappelt” — sie pendelte sich Runde um Runde in die richtige Antwort ein. Das war erwünscht: der Lösungsweg selbst. Hier ist das Zappeln das genaue Gegenteil — es ist ein Verfahrensfehler, der sich aufschaukelt und von der Wahrheit wegläuft. Beim Einpendeln werden die Ausschläge kleiner, bei der Explosion größer. Dieselbe Zickzack-Optik, entgegengesetzte Bedeutung.

9.4.1 Die scharfe Grenze und die Viertel-Regel

Warum kippt es genau zwischen 1800 und 3000 Sekunden? Der explizite Schritt korrigiert jeden Knoten um einen Bruchteil seines Rückstands zu den Nachbarn. Ist dieser Bruchteil zu groß, schießt die Korrektur über das Ziel hinaus — und weil im nächsten Schritt dasselbe mit umgekehrtem Vorzeichen passiert, wächst der Überschuss. Die Grenze, ab der das losgeht, lässt sich als Faustregel hinschreiben:

\[ \Delta t_{\max} \;=\; \frac{C_{\text{innen}}}{2\,k} \;=\; \frac{9027{,}5}{2\cdot 2} \;=\; \frac{9027{,}5}{4} \;\approx\; 2257\ \mathrm{s} . \]

Das sind rund 38 Minuten. Unsere 1800 s (eine halbe Stunde) liegen darunter (stabil), die 3000 s (fünfzig Minuten) darüber (Explosion) — genau wie beobachtet. Dieselbe Grenze lässt sich auch mit den Materialgrößen und der Maschenweite schreiben, und in dieser Form zeigt sie das Entscheidende:

\[ \Delta t_{\max} \;=\; \frac{h^2\,\rho\,c}{2\,\lambda} \;=\; \frac{0{,}25^2 \cdot 3\,611\,000}{2\cdot 50} \;\approx\; 2257\ \mathrm{s} . \]

Das \(h^2\) ist der Haken. Die Maschenweite steht im Quadrat. Halbiert man sie — verdoppelt also die Knotenzahl, um ein feineres Bild zu bekommen —, dann viertelt sich der erlaubte Zeitschritt. Rechnen wir es für \(h = 0{,}125\ \mathrm{m}\) (9 statt 5 Knoten) nach:

\[ \Delta t_{\max} \;=\; \frac{0{,}125^2 \cdot 3\,611\,000}{2\cdot 50} \;\approx\; 564\ \mathrm{s} \qquad(\text{ein Viertel von }2257\ \mathrm{s}). \]

Code
import numpy as np
import matplotlib.pyplot as plt

rhoc, lam = 3_611_000.0, 50.0
h = np.linspace(0.05, 0.5, 200)
dt_max = h**2 * rhoc / (2 * lam)

fig, ax = plt.subplots(figsize=(8.0, 4.6))
ax.fill_between(h, 0, dt_max, color="#d9f0d3", label="stabil")
ax.fill_between(h, dt_max, 9500, color="#fdecea", label="explodiert")
ax.plot(h, dt_max, "-", color="black", lw=2.0)

# Die Faustregelpunkte fuer h = 0,25 und h = 0,125.
ax.plot(0.25, 2257, "s", color="black", ms=8)
ax.annotate("h = 0,25 m → 2257 s", xy=(0.25, 2257), xytext=(0.30, 3400),
            fontsize=9, arrowprops=dict(arrowstyle="->"))
ax.plot(0.125, 564, "s", color="black", ms=8)
ax.annotate("h = 0,125 m → 564 s", xy=(0.125, 564), xytext=(0.13, 2200),
            fontsize=9, arrowprops=dict(arrowstyle="->"))

# Die beiden Experimente bei h = 0,25.
ax.plot(0.25, 1800, "o", color="tab:green", ms=11, markeredgecolor="black")
ax.text(0.265, 1500, "Δt = 1800 s\n(stabil)", fontsize=9, color="tab:green")
ax.plot(0.25, 3000, "o", color="tab:red", ms=11, markeredgecolor="black")
ax.text(0.265, 3050, "Δt = 3000 s\n(Explosion)", fontsize=9, color="tab:red")

ax.set_xlim(0.05, 0.5)
ax.set_ylim(0, 9500)
ax.set_xlabel("Maschenweite h (m)")
ax.set_ylabel("erlaubter Zeitschritt Δt (s)")
ax.set_title("Stabilitätskarte des expliziten Euler")
ax.legend(loc="upper left")
plt.tight_layout()
plt.show()
Abbildung 9.6: Die Stabilitätskarte des expliziten Euler: erlaubter Zeitschritt gegen Maschenweite. Unter der Kurve (grün) ist der Schritt stabil, darüber (rot) explodiert er. Die Kurve fällt mit dem Quadrat der Maschenweite. Eingezeichnet: die beiden Experimentpunkte bei h = 0,25 m — 1800 s liegt sicher im Grünen, 3000 s im Roten.

Was man hier sieht: eine Landkarte mit einer grünen und einer roten Zone, getrennt durch die Faustregel-Kurve. Wer ein feineres Netz will, rutscht auf der Kurve nach links — und der erlaubte Zeitschritt stürzt überproportional ab. Ein doppelt so feines Netz braucht nicht doppelt, sondern viermal so viele Zeitschritte, und jeder ist auch noch teurer. Feine explizite Rechnungen werden so schnell unbezahlbar. Das ist der Grund, aus dem es das nächste Verfahren gibt.

9.5 Die implizite Rettung

Der explizite Schritt fragt: „Wie schnell ändern sich die Dinge jetzt?” und schreibt diese Rate für den ganzen Schritt fest. Das ist der Fehler: Am Anfang des Schritts ist die Rate am größten, und wenn der Schritt zu lang ist, überschießt sie. Der implizite Euler-Schritt stellt die Frage anders: „Welche neuen Werte passen so zusammen, dass die Änderungs- rate am Ende des Schritts stimmt?”

Das klingt nach einer Henne-Ei-Frage — die neue Rate hängt von den neuen Werten ab, die wir doch erst suchen. Und genau das ist der Punkt: Man kann die Knoten nicht mehr einzeln nacheinander ausrechnen, sondern nur noch alle gleichzeitig. Für jeden inneren Knoten \(i\) lautet die Forderung

\[ \Bigl(\frac{C_i}{\Delta t} + 2k\Bigr) T_i^{\text{neu}} \;-\; k\,T_{i-1}^{\text{neu}} \;-\; k\,T_{i+1}^{\text{neu}} \;=\; \frac{C_i}{\Delta t}\,T_i^{\text{alt}} . \]

Jede Zeile koppelt einen Knoten an seine zwei Nachbarn — das ist ein tridiagonales Gleichungssystem, dieselbe Bandstruktur wie beim stationären Träger in Kapitel 7. Und dafür haben wir längst das perfekte Werkzeug: den Thomas-Algorithmus aus Kapitel 4. Ein impliziter Zeitschritt ist also ein Thomas-Lauf. Er kostet mehr als die simple Division des expliziten Schritts — aber er hat eine Eigenschaft, die alles aufwiegt: Er ist bei jedem Zeitschritt stabil. Kein \(\Delta t\) ist zu groß; die Explosion kann nicht passieren.

WarnungVermutung aufgelöst

Die Vermutung vom Kapitelanfang — „der Zeitschritt ist nur eine Frage der Genauigkeit” — ist beim expliziten Verfahren schlicht falsch: Jenseits der scharfen Grenze explodiert es. Beim impliziten Verfahren wird sie dagegen fast wahr: Stabil ist es immer, und ein größerer Schritt macht das Ergebnis nur ungenauer, nicht kaputt. „Fast”, weil ungenau hier etwas Bestimmtes heißt — die Front verschmiert. Das schauen wir uns als Letztes an.

9.5.1 Stabil heißt nicht genau

Man könnte nun übermütig werden: Wenn implizit immer stabil ist, nehme ich eben einen riesigen Zeitschritt und bin in einem Sprung fertig. Stabil bleibt es — aber genau bleibt es nicht. Grobe Schritte verschmieren die Front: Die scharfe Kante zwischen warm und kalt, die durch den Träger wandert, wird weichgezeichnet.

Code
import numpy as np
import matplotlib.pyplot as plt

lam, A, L = 50.0, 0.01, 1.0
rho, c = 7850.0, 460.0
knoten = 21
h = L / (knoten - 1)
k = lam * A / h
C = rho * c * A * h
orte = [i * h for i in range(knoten)]

def implizit(dt, endzeit):
    haupt = (C / dt + 2 * k) * np.ones(knoten - 2)
    neben = -k * np.ones(knoten - 3)
    M = np.diag(haupt) + np.diag(neben, 1) + np.diag(neben, -1)
    T = np.full(knoten, 20.0); T[-1] = 300.0
    for s in range(int(endzeit / dt)):
        rs = (C / dt) * T[1:-1].copy()
        rs[0] += k * 20.0
        rs[-1] += k * 300.0
        T[1:-1] = np.linalg.solve(M, rs)
    return T

endzeit = 3600.0
referenz = implizit(60.0, endzeit)
mittel = implizit(600.0, endzeit)
grob = implizit(3600.0, endzeit)   # ein einziger Riesenschritt

fig, ax = plt.subplots(figsize=(8.5, 4.4))
ax.plot(orte, referenz, "-", color="black", lw=2.2, label="feine Referenz")
ax.plot(orte, mittel, "o-", color="tab:blue", ms=5,
        label="implizit, Δt = 600 s")
ax.plot(orte, grob, "s--", color="tab:purple", ms=6,
        label="implizit, Δt = 3600 s (1 Schritt)")
ax.set_xlabel("Ort entlang des Trägers (m)")
ax.set_ylabel("Temperatur (°C)")
ax.set_title("Nach 3600 s: gleich stabil, ungleich genau")
ax.legend()
plt.tight_layout()
plt.show()
Abbildung 9.7: Genauigkeit ist nicht Stabilität. Alle drei Kurven sind nach einer Stunde (3600 s) auf demselben 21-Knoten-Träger gerechnet und alle drei stabil. Die feine Referenz (schwarz, winzige Schritte) zeigt die echte, noch recht scharfe Front. Der implizite Lauf mit mittlerem Δt (blau) trifft sie gut; der eine Riesenschritt (lila) hat die Front verschmiert — dieselbe Endzeit, dieselbe Stabilität, weniger Genauigkeit.

Was man hier sieht: drei stabile Rechnungen, ein und dieselbe Endzeit. Die feine Referenz zeigt die echte Front — rechts schon deutlich warm, links noch kalt, mit einem klaren Übergang. Der mittlere Zeitschritt (blau) liegt fast auf ihr. Der eine Riesenschritt (lila) dagegen hat die Front verschmiert: Er unterschätzt die Wärme rechts und schmiert sie zu weit nach links. Kein Absturz, kein Zappeln — nur eben ungenau. Genau davor warnt die aufgelöste Vermutung: Stabilität ist geschenkt, Genauigkeit nicht.

9.5.2 Explizit oder implizit — was nehmen?

Code
import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10.0, 3.4))
ax.axis("off")
zeilen = [
    ["", "explizit", "implizit"],
    ["Kosten je Schritt", "eine Division je Knoten",
     "ein Gleichungssystem (Thomas)"],
    ["Zeitschritt", "an die Maschenweite gefesselt (h²)", "frei wählbar"],
    ["Stabilität", "nur unter Δt_max", "immer stabil"],
    ["Gefahr bei groß Δt", "Explosion", "nur Front verschmiert"],
    ["gut, wenn …", "Netz grob, Schritt klein gewollt",
     "Netz fein oder große Schritte gewollt"],
]
tab = ax.table(cellText=zeilen, loc="center", cellLoc="left")
tab.auto_set_font_size(False)
tab.set_fontsize(10)
tab.scale(1.0, 1.9)
for c in range(3):
    tab[0, c].set_facecolor("#d9e6f2")
    tab[0, c].set_text_props(weight="bold")
for r in range(len(zeilen)):
    tab[r, 0].set_text_props(weight="bold")
plt.tight_layout()
plt.show()
Abbildung 9.8: Die Entscheidungstabelle: explizit gegen implizit. Der explizite Schritt ist billig, aber sein Zeitschritt ist an die Netzfeinheit gefesselt. Der implizite Schritt kostet pro Schritt ein Gleichungssystem, darf dafür aber so große Schritte machen, wie es die gewünschte Genauigkeit erlaubt.

Was man hier sieht: kein klarer Sieger, sondern ein Tausch. Wer ein grobes Netz hat und ohnehin kleine Schritte will, fährt mit dem billigen expliziten Verfahren gut. Wer ein feines Netz braucht (und darum an die \(h^2\)-Fessel gerät) oder große Schritte machen möchte, nimmt das implizite: Sein Mehraufwand je Schritt zahlt sich aus, weil er unvergleichlich weniger Schritte braucht. Genau deshalb rechnen die meisten professionellen Programme transiente Wärme implizit — die Stabilität ist ihnen den Preis wert.

9.6 Zwanzig Stunden — und der Brandschutz rechnet in Minuten

Eine Zahl fehlt noch, und sie erklärt, warum in jeder Animation dieses Kapitels „Zeitraffer” steht. Wie lange dauert das Aufheizen überhaupt? Die Antwort steckt in einer einzigen Kombination der Trägergrößen, der Zeitkonstante:

\[ \tau \;=\; \frac{L^2\,\rho\,c}{\lambda} \;=\; \frac{1{,}0^2 \cdot 3\,611\,000}{50} \;=\; 72\,220\ \mathrm{s} \;\approx\; 20\ \mathrm{Stunden} . \]

HinweisEhrlich gerechnet: die Einheiten

\(L^2\) hat die Einheit \(\mathrm{m^2}\), \(\rho c\) ist \(\mathrm{J/(m^3\,K)}\), \(\lambda\) ist \(\mathrm{W/(m\,K)} = \mathrm{J/(s\,m\,K)}\). Alles zusammengesetzt kürzen sich Meter, Joule und Kelvin weg, und übrig bleiben Sekunden — wie es sich für eine Zeit gehört. Das ist keine erfundene Formel, sondern die natürliche Zeitskala, die aus der Wärmeleitungsgleichung selbst fällt.

Rund zwanzig Stunden also, bis unser Stahlträger über seine ganze Länge eingeschwungen ist. In Echtzeit wäre das eine sehr langweilige Animation. Darum lief der Film am Kapitelanfang im Zeitraffer — zehn Stunden in wenigen Sekunden, mit der eingeblendeten Uhr, die die echten Stunden zählt, ganz wie die Zeitraffer-Uhr aus Kapitel 1. Ehrlich beschriftet ist der Zeitraffer kein Trick, sondern eine Notwendigkeit.

Und genau hier steckt die eigentliche Pointe für den Brandschutz — sie verlangt einen eigenen Kasten.

WichtigR30, R60, R90 — Brandschutz rechnet in Minuten

Bauteile werden nach ihrer Feuerwiderstandsdauer eingeteilt: R30, R60, R90 heißt, das Bauteil hält im Normbrand 30, 60 oder 90 Minuten, bevor es seine Tragfähigkeit verliert (CEN (Europäisches Komitee für Normung) 2005). Minuten — nicht Stunden. Wie weit ist die Wärme in dieser Zeit überhaupt in den Träger vorgedrungen? Das schätzt die Eindringtiefe \(\sqrt{a\,t}\) mit der Temperaturleitfähigkeit \(a = \lambda/(\rho c) = 1{,}385\cdot10^{-5}\ \mathrm{m^2/s}\):

\[ \sqrt{a\cdot 30\ \mathrm{min}} = \sqrt{1{,}385\cdot10^{-5}\cdot 1800} \;\approx\; 0{,}16\ \mathrm{m} \;=\; 16\ \mathrm{cm} . \]

Nach einer halben Stunde ist die Hitze also erst rund 16 cm weit gekrochen, nach 90 Minuten (R90) rund 27 cm — der ganze Träger ist einen Meter lang. Das kühle Auflager am fernen Ende bleibt in der gesamten brandschutzrelevanten Zeit kalt. Die Zeitkonstante \(\tau \approx 20\ \mathrm{h}\) sagt, wann der ganze Träger durchgewärmt wäre; der Brandschutz aber entscheidet sich in den ersten Minuten, nahe dem Brandherd. Beide Uhren gehören zusammen — und beide muss man kennen.

HinweisEin Vorgriff: warum später 1D genügt

Dieselbe Zeitkonstante lässt sich auch über die Höhe des Trägers bilden. Er ist nur \(H = 0{,}1\ \mathrm{m}\) hoch, und \(H\) geht im Quadrat ein:

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

Quer über die Höhe gleicht sich ein Temperaturunterschied also in etwa 12 Minuten aus — rund hundertmal schneller als über die Länge (20 Stunden). Genau deshalb dürfen wir die Wärme im Träger überhaupt als 1D-Problem entlang der Länge rechnen: quer ist der Träger praktisch immer schon ausgeglichen. Diesen Gedanken nehmen wir in Kapitel 10 wieder auf, wenn wir in die Fläche gehen.

9.7 Die interaktive Einheit: die Zeitmaschine

Jetzt drehst du selbst an der Uhr. Die Zeitmaschine rechnet den transienten Träger mit beiden Verfahren; du stellst Zeitschritt, Knotenzahl und Verfahren ein und siehst zweierlei: das Temperaturprofil am Ende und die Temperatur des Mittenknotens über der ganzen Zeit. Der Thomas-Löser aus Kapitel 4 ist unten hineinkopiert, weil der Browser nicht aus den Kapitel-Ordnern importieren kann.

WichtigVorhersage-Punkt

Bevor du ausführst: Die Zeitmaschine startet mit 5 Knoten, explizit und \(\Delta t = 1800\ \mathrm{s}\) — sauber unter der Grenze von 2257 s, also stabil. Ändere jetzt in der ersten Zeile ZEITSCHRITT_S = 1800.0 auf 3000.0 und führe erneut aus. Was zeigt die Kurve des Mittenknotens dann — läuft sie glatt auf ihren Endwert, oder tut sie etwas anderes? Leg dich fest, bevor du es probierst.

Abbildung 9.9: Vorgerenderte Fassung der Zeitmaschine mit der Startbelegung (5 Knoten, explizit, Δt = 1800 s): links das eingeschwungene Profil, rechts der Mittenknoten, der glatt auf seine 160 °C zuläuft. Im Browser wird dieses Bild durch dein eigenes Ergebnis ersetzt.

Was man hier sieht: links das vertraute Endprofil, die Gerade von 20 auf 300 °C, mit den rund 160 °C in der Mitte; rechts der Weg dorthin — der Mittenknoten steigt erst zögernd (die Wärme ist noch nicht da), dann zügig, dann immer flacher, bis er sich bei 160 °C anlegt. Stell in der ersten Zeile ZEITSCHRITT_S = 3000.0 ein, und diese glatte Kurve wird zur zackigen Explosion. Schalte danach auf VERFAHREN = "implizit" um — und selbst mit \(\Delta t = 3000\ \mathrm{s}\) ist wieder Ruhe.

9.8 Das Kapitel-Programm

Das vollständige, eigenständige Programm zu diesem Kapitel liegt in programme/kap09/kap09_transiente_fem.py. Es baut die gelumpten Kapazitäten, bietet ein_expliziter_schritt und ein_impliziter_schritt (letzterer importiert den Thomas-Löser aus programme/gemeinsam/loeser.py) und rechnet beim Aufruf die Handrechnung, die Explosion und die Frontverschmierung nach. Die im Kapitel gezeigten Zahlen stammen alle aus diesem Programm.

TippMerkkasten
  • Der Speicherterm \(\rho c\,\partial T/\partial t\) ist die Bilanz aus Kapitel 6, bei der die Differenz zwischen Zu- und Abfluss nicht mehr null ist, sondern das Stück aufheizt. \(\rho c = 3{,}61\ \mathrm{MJ/(m^3\,K)}\) ist der Preis der Trägheit.
  • Der Speicherterm wird zur Massenmatrix. Lumping macht sie diagonal: Jeder Knoten bekommt seinen Eimer — innen 9027,5 J/K, am Rand 4513,75 J/K.
  • Expliziter Euler: neuer Wert = alter Wert + Zeitschritt · Änderungs- rate. Billig (eine Division je Knoten), aber nur unter \(\Delta t_{\max} = h^2\rho c/(2\lambda) \approx 2257\ \mathrm{s}\) (rund 38 min) stabil — darüber explodiert die Rechnung.
  • Halbe Maschenweite → ein Viertel des Zeitschritts (das \(h^2\)).
  • Impliziter Euler: alle neuen Werte zugleich aus einem tridiagonalen System (Thomas). Immer stabil, aber grobe Schritte verschmieren die Front. Stabilität ≠ Genauigkeit.

Rückblick auf die Landkarte

Mit diesem Kapitel endet Teil III — und der Baustein Physik unserer Landkarte ist damit vollständig. Wir sind von der Steigung und der Fläche (Kapitel 5) über die Bilanz (Kapitel 6) und die schwache Form (Kapitel 7) bis zu den Randbedingungen (Kapitel 8) gekommen, und jetzt tickt auch die Uhr (dieses Kapitel). Der Träger kann stationär und transient gerechnet werden.

Code
import matplotlib.pyplot as plt

bausteine = [
    ("Geometrie\n& Netz", "Kap. 10–11", False),
    ("Physik", "Kap. 5, 6, 9", True),
    ("Rand-\nbedingungen", "Kap. 8", True),
    ("Gleichungs-\nsystem", "Kap. 3, 7", True),
    ("Löser", "Kap. 4", True),
]
farben = ["#cfe8f3", "#ffe0b3", "#ffd6d6", "#d9f0d3", "#e6dcf0"]

fig, ax = plt.subplots(figsize=(9.0, 3.0))
breite, luecke = 1.5, 0.4
for k, ((name, kapitel, fertig), farbe) in enumerate(zip(bausteine, farben)):
    x0 = k * (breite + luecke)
    rand = "black" if fertig else "#bbbbbb"
    breite_rand = 2.6 if fertig else 1.0
    alpha = 1.0 if fertig else 0.45
    ax.add_patch(plt.Rectangle((x0, 0), breite, 1.0, facecolor=farbe,
                               edgecolor=rand, lw=breite_rand, alpha=alpha))
    ax.text(x0 + breite / 2, 0.62, name, ha="center", va="center",
            fontsize=10, weight="bold")
    ax.text(x0 + breite / 2, 0.22, kapitel, ha="center", va="center",
            fontsize=9)
    if fertig:
        ax.text(x0 + breite / 2, 0.90, "✓", ha="center", va="center",
                fontsize=13, color="tab:green", weight="bold")
    if k < len(bausteine) - 1:
        ax.annotate("", xy=(x0 + breite + luecke, 0.5),
                    xytext=(x0 + breite, 0.5),
                    arrowprops=dict(arrowstyle="-|>", color="black", lw=1.3))

ax.text((len(bausteine) * (breite + luecke) - luecke) / 2, 1.45,
        "Teil III abgeschlossen — als Nächstes: in die Fläche (Teil IV)",
        ha="center", va="center", fontsize=9, style="italic")
ax.set_xlim(-0.2, len(bausteine) * (breite + luecke) - luecke + 0.2)
ax.set_ylim(-0.2, 1.7)
ax.axis("off")
plt.tight_layout()
plt.show()
Abbildung 9.10: Die Landkarte des Buches am Ende von Teil III: Der Baustein „Physik” ist jetzt komplett ausgefüllt (kräftig hervorgehoben) — stationär und transient. Der Löser (Kap. 4) steht ebenfalls. Was noch fehlt, ist die Fläche: Geometrie und Netz jenseits der 1D-Kette (Teil IV).

Was man hier sieht: dieselbe Landkarte wie in Kapitel 1, aber vier von fünf Kästen tragen jetzt ihr Häkchen. Nur der erste — Geometrie und Netz — ist noch blass: Bisher war unser „Netz” eine simple Kette aus fünf Knoten auf einer Linie. Teil IV bricht in die Fläche auf, mit Dreiecken und Rechtecken, und füllt auch diesen letzten Kasten.

Roter Faden

Zurück: Der Speicherterm ist die Bilanz aus Kapitel 6, diesmal mit der Differenz, die aufheizt. Die Massenmatrix wird nach derselben Halbe-Halbe-Logik gebaut wie die Wärmequelle in Kapitel 7. Der implizite Schritt ist ein tridiagonales System und damit ein Fall für den Thomas-Algorithmus aus Kapitel 4. Die festen Randtemperaturen behandeln wir wie in Kapitel 8. Und die explizite Zeitschleife aus der Blackbox von Kapitel 1 ist jetzt kein Geheimnis mehr.

Vor: Lumping und Zeitschritte laufen in 2D unverändert weiter (Kapitel 10 und Kapitel 15). Die Explosion ist der Grund, aus dem das große Finale (Kapitel 15) implizit rechnet. Und die Frontverschmierung wird in Kapitel 12 zu einem Prüffall: Woher weiß man, ob ein „schön aussehendes” Ergebnis nicht bloß eine grob verschmierte Näherung ist?

Übungen

Ü 9.1 (Verstehen). Führe die Handrechnung aus Abbildung 9.4 um einen vierten Schritt fort (von „nach 1800 s” auf „nach 2400 s”), mit denselben Zahlen \(\Delta t = 600\ \mathrm{s}\), \(k = 2\ \mathrm{W/K}\), \(C = 9027{,}5\ \mathrm{J/K}\). Rechne die drei inneren Knoten von Hand und prüfe mit der Zeitmaschine (5 Knoten, explizit, \(\Delta t = 600\), 4 Schritte).

Vorfaktor wie gehabt \(\frac{600}{9027{,}5}\cdot 2 \approx 0{,}1329\). Aus dem Zustand nach 1800 s, \([20;\ 20{,}658;\ 32{,}212;\ 105{,}262;\ 300]\):

  • Knoten 2: \(20{,}658 + 0{,}1329\cdot(20 - 2\cdot 20{,}658 + 32{,}212) = 22{,}106\ °\mathrm{C}\),
  • Knoten 3: \(32{,}212 + 0{,}1329\cdot(20{,}658 - 2\cdot 32{,}212 + 105{,}262) = 40{,}386\ °\mathrm{C}\),
  • Knoten 4: \(105{,}262 + 0{,}1329\cdot(32{,}212 - 2\cdot 105{,}262 + 300) = 121{,}438\ °\mathrm{C}\).

Ergebnis nach 2400 s: \([20;\ 22{,}11;\ 40{,}39;\ 121{,}44;\ 300]\). Die Wärme kriecht weiter nach links; Knoten 4 hat gerade die 120 °C überschritten.

Ü 9.2 (Verändern). Verdopple die Zahl der Elemente (5 → 9 Knoten, \(h\) von 0,25 auf 0,125 m) und kreise mit der Zeitmaschine den größten noch stabilen Zeitschritt experimentell ein: Erhöhe \(\Delta t\) Schritt für Schritt, bis die Kurve des Mittenknotens explodiert. Vergleiche mit der Faustregel. Stimmt die Viertel-Regel?

Die Faustregel liefert für 9 Knoten \(\Delta t_{\max} = h^2\rho c/(2\lambda) = 0{,}125^2\cdot 3\,611\,000/100 \approx 564\ \mathrm{s}\) — ein Viertel der 2257 s des 5-Knoten-Trägers, weil \(h\) im Quadrat steht. Experimentell bleibt der Lauf bis knapp über diese Grenze ruhig und kippt kurz darüber. Das Lösungsskript loesungen/kap09_ue2.py kreist beide Grenzen ein und bestätigt das Verhältnis 4 : 1.

Ü 9.3 (Übertragen). Tausche das Material von Baustahl auf Aluminium (\(\lambda = 200\ \mathrm{W/(m\,K)}\), \(\rho = 2700\ \mathrm{kg/m^3}\), \(c = 900\ \mathrm{J/(kg\,K)}\)). Rechne zuerst von Hand die Zeitkonstante \(\tau = L^2\rho c/\lambda\) und die Faustregel \(\Delta t_{\max}\) aus, dann simuliere. Wird ein Aluminiumträger schneller oder langsamer warm — und darf der explizite Schritt größer oder kleiner sein als bei Stahl?

\(\rho c = 2700\cdot 900 = 2{,}43\ \mathrm{MJ/(m^3\,K)}\), also \(\tau = 1\cdot 2\,430\,000/200 = 12\,150\ \mathrm{s} \approx 3{,}4\ \mathrm{Stunden}\) — fast sechsmal flinker als Stahl (Aluminium leitet die Wärme viel besser). Die Faustregel: \(\Delta t_{\max} = 0{,}25^2\cdot 2\,430\,000/(2\cdot 200) \approx 380\ \mathrm{s}\), also fast sechsmal kleiner als bei Stahl — das große \(\lambda\) steht im Nenner. Beides folgt aus derselben Ursache: Aluminium leitet gut, also flitzt die Wärme schnell und muss in kleineren Zeitsprüngen verfolgt werden. Das Lösungsskript loesungen/kap09_ue3.py rechnet und simuliert das.

Das Kleingedruckte

Woher die konsistente Massenmatrix kommt. Die \(6\) und die \(2,1\) im Ergebnis-Kasten sind die exakten Flächen unter den Produkten zweier Hütchen. Über ein Element der Länge \(h\) ist \(\int N_i N_i\,\mathrm{d}x = h/3\) (die Fläche unter einem quadratischen Hütchenquadrat) und \(\int N_i N_j\,\mathrm{d}x = h/6\) für Nachbarn. Mal \(\rho c A\) ergibt das \(\frac{\rho c A h}{6}\left[\begin{smallmatrix} 2 & 1 \\ 1 & 2\end{smallmatrix}\right]\). Das Lumping addiert je Zeile zusammen und schreibt die Summe auf die Diagonale — massenerhaltend, wie die Eimer zeigen. Wer die Flächen selbst nachrechnen will, findet das Werkzeug in Kapitel 5.

Was dieses Kapitel bewusst weglässt. Es gibt ein Verfahren genau zwischen explizit und implizit — das Crank-Nicolson-Verfahren, das die Änderungsrate zur Hälfte am Anfang und zur Hälfte am Ende ansetzt und dadurch genauer wird. Es ist der übliche Kompromiss in der Praxis; wir erwähnen es nur. Ebenso haben wir die Stabilitätsgrenze als Faustregel benutzt, nicht als exakten Beweis. Für die kleine 5-Knoten-Kette liegt die echte Grenze etwas höher (rund 2640 s, also gut 44 Minuten) als die 2257 s (rund 38 Minuten) der Faustregel — die Faustregel ist bewusst auf der sicheren Seite. Wer es genau wissen will, braucht die sogenannte von-Neumann-Stabilitätsanalyse; sie steht nicht auf dem Weg dieses Buches. Auch die feine Frage, ob die konsistente Massenmatrix bei sehr kleinen Zeitschritten leicht überschwingt, lassen wir beim Eimerbild bewenden: Gelumpt ist robust, und robust genügt uns.

CEN (Europäisches Komitee für Normung). 2005. Eurocode 3: Bemessung und Konstruktion von Stahlbauten – Teil 1-2: Allgemeine Regeln – Tragwerksbemessung für den Brandfall. EN 1993-1-2. CEN.
Euler, Leonhard. 1768. Institutiones calculi integralis. Bd. 1. Impensis Academiae Imperialis Scientiarum.