Jede dynamische Gebäudesimulation steht und fällt mit ihrem Rechenkern — dem Stück Software, das aus Geometrie, Bauteilen, Wetter und Nutzung einen physikalisch belastbaren Temperatur- und Lastverlauf macht. Bei VICUS Buildings ist dieser Kern NANDRAD, ein an der TU Dresden in über 20 Jahren Forschung gereifter Simulationslöser. Dieser Artikel öffnet die Motorhaube: Er zeigt Schritt für Schritt, was der NandradSolver rechnet und wie er es tut — von der physikalischen Bilanzgleichung über die Zeitintegration bis zum linearen Gleichungssystem im Inneren.

Sie müssen dafür keine Numerik studiert haben. Die zentrale Idee lässt sich in einem Satz sagen und mit den interaktiven Demos weiter unten selbst „durchkurbeln”.

Die Grundidee in einem Satz

Der Solver beschreibt das Gebäude (und gegebenenfalls ein Rohrnetz) als ein großes System gekoppelter gewöhnlicher Differentialgleichungen und integriert dieses System Zeitschritt für Zeitschritt durch die Zeit. Das Herzstück ist immer dieselbe Frage:

„Wie schnell ändert sich die gespeicherte Energie in jedem Raum und in jeder Wandschicht?”

Kennt man diese Änderungsrate, kann man einen kleinen Schritt in die Zukunft gehen, die neuen Temperaturen berechnen und das Spiel von vorne beginnen. In Vektorschreibweise — und genau so „sieht” es der Zeitintegrator — lautet das:

$\frac{d\mathbf{y}}{dt} = \texttt{ydot}(t,\,\mathbf{y})$

Dabei ist $\mathbf{y}$ der Zustandsvektor (eine lange Liste aller Erhaltungsgrößen), $\texttt{ydot}$ deren zeitliche Änderung und $t$ die Zeit.

Das Schema und drei Live-Demos

Die folgende interaktive Visualisierung fasst den gesamten Ablauf zusammen — und macht ihn anfassbar. Oben das Funktionsschema mit den drei geschachtelten Schleifen (fahren Sie über die Boxen, um die zugehörigen Quelldateien zu sehen; schalten Sie zwischen implizit und explizit um). Darunter drei Live-Demos, in denen Sie selbst an Zeitschritt, Netzfeinheit und Wärmetauscher-Parametern drehen und die Stabilitätsgrenze finden können.

Der Rest des Artikels erklärt die einzelnen Bausteine im Detail. Sie können jederzeit nach oben scrollen und das Gelernte direkt ausprobieren.

Energie statt Temperatur: die Erhaltungsgröße

Eine zentrale Designentscheidung in NANDRAD: Der Solver bilanziert Energie (genauer: die gespeicherte innere Energie), nicht direkt die Temperatur. Der Grund ist physikalisch sauber — Energie ist eine Erhaltungsgröße. Man kann alle Energieströme, die in einen Raum hinein- und hinausfließen, eindeutig aufsummieren. Die Temperatur ist dann eine abgeleitete Größe, die man aus der Energie und der Wärmekapazität zurückrechnet.

Die allgemeine Form jeder Bilanzgleichung lautet:

$\frac{d(\text{gespeicherte Energie})}{dt} = \sum \text{Zuflüsse} - \sum \text{Abflüsse}$

Was steckt im Zustandsvektor y?

y ist schlicht alles aneinandergehängt, was sich über die Zeit ändert und Energie speichert:

Eine einzige Wand wird also intern in viele dünne Schichten zerlegt — jede Schicht ist ein eigener Eintrag in y. Bei einem größeren Gebäude kommen so schnell mehrere tausend Unbekannte zusammen. Genau diese Zerlegung „einer Wand in viele gekoppelte Gleichungen” können Sie in der zweiten Demo (Wärmeleitung durch eine Wand) am Schieberegler N selbst erzeugen.

Woher kommt ydot? Die Bilanz, nicht geraten

Eine berechtigte Frage: $dy/dt = \texttt{ydot}(t, y)$ steht da — aber wie kommt man auf die rechte Seite? Antwort: ydot rät man nicht, man leitet es aus einem Erhaltungssatz her. Das Rezept ist immer dasselbe:

1. Wähle die Erhaltungsgröße y (Energie, Masse, …) — nicht die Temperatur.
2. Schreibe die Bilanz: d(Gespeichertes)/dt = Σ Zuflüsse − Σ Abflüsse.
3. Drücke jeden Fluss durch den Zustand aus (physikalische Gesetze).
4. Löse nach dy/dt auf — das ist ydot.

Beispiel: ein Raum kühlt aus

Gespeichert wird die innere Energie $U = C\,T$ mit der Wärmekapazität $C$. Sie ändert sich nur durch den Wärmestrom $Q$ über die Hülle, und dieser ist proportional zur Temperaturdifferenz:

$\frac{dU}{dt} = Q, \qquad Q = -\frac{T - T_\text{außen}}{R}$

Mit $\frac{dU}{dt} = C\,\frac{dT}{dt}$ folgt durch Auflösen nach der Rate:

$\frac{dT}{dt} = -\frac{T - T_\text{außen}}{R\,C} = -\frac{T - T_\text{außen}}{\tau}, \qquad \tau = R\,C$

Das ydot ist damit hergeleitet, nicht erfunden — Energieerhaltung plus ein Flussgesetz. Genau diese eine Gleichung läuft in der ersten Demo („ein Raum kühlt aus”): Sie sehen, wie der Integrator sich an der Tangente entlang durch die Zeit hangelt und mit der analytischen Lösung verglichen wird.

Und bei einem echten Gebäude?

Die entscheidende Einsicht: In NANDRAD gibt es kein einziges geschlossenes ydot-Formelblatt. Bei tausenden Zuständen wäre eine Hand-Formel unmöglich. Stattdessen ist die Bilanz dieselbe für jeden Zustand, aber die Summe der Flüsse wird zur Laufzeit zusammengesetzt — jedes Teilmodell wirft seinen Beitrag in den Topf: Wärmeleitung zwischen Schichten, solare Einstrahlung durch Fenster, Lüftung und Infiltration, interne Lasten, Heiz- und Kühlleistung, Rohrströmung. ydot ist deshalb eine Funktion, die man nicht hinschreibt, sondern auswertet.

Die drei Auswertungsstufen: Head → Mitte → Tail

Der Integrator übergibt dem Modell ein y und fragt nach dem zugehörigen ydot. Intern läuft das in drei Stufen ab — genau das Bilanz-Rezept, in Code gegossen:

Im Schema oben ist dieser Datenfluss als dreistufige Pipeline innerhalb der Newton-Schleife dargestellt.

Reihenfolge der Modelle: der Abhängigkeitsgraph

Es gibt viele kleine Teilmodelle, und sie hängen voneinander ab: Das Fenster-Modell braucht die Außentemperatur vom Klima-Modell, die Raumbilanz die Heizleistung vom HVAC-Modell. In welcher Reihenfolge muss man rechnen?

NANDRAD löst das mit einem Abhängigkeitsgraphen. Jedes Modell meldet an, welche Größen es als Eingabe braucht und welche es als Ergebnis liefert. Daraus wird automatisch eine korrekte Auswertungsreihenfolge berechnet (topologische Sortierung). Zwei Spezialfälle sind bemerkenswert:

• Parallelisierung: Modelle, die voneinander unabhängig sind, liegen in derselben „Ebene” des Graphen und werden parallel ausgewertet (z. B. viele Wände gleichzeitig).
• Zyklen: Hängen Modelle gegenseitig voneinander ab (etwa ein Regler, der auf eine Größe reagiert, die er selbst beeinflusst), werden sie zu einer Gruppe zusammengefasst und iterativ gemeinsam gelöst.

Zeitintegration: einen Schritt in die Zukunft gehen

Jetzt haben wir eine Funktion, die zu jedem y das ydot liefert. Die Aufgabe der Zeitintegration ist, daraus den zeitlichen Verlauf zu rekonstruieren — aus der Steigung den Weg. In einer äußeren Schleife geht der Solver Schritt für Schritt von t nach t+dt, schreibt jeden akzeptierten Schritt fort und schreibt fällige Ergebnisse heraus.

Wichtig: Die Schrittweite dt ist nicht konstant. Der Integrator wählt sie selbst — mal große Schritte (nachts, wenn sich wenig ändert), mal sehr kleine (wenn die Heizung anspringt). Die Ausgabe (etwa stündliche Werte) wird davon entkoppelt: Zwischen den tatsächlichen Rechenschritten wird sauber auf die gewünschten Ausgabezeitpunkte interpoliert.

Welcher Integrator? CVODE als Standard

NANDRAD bringt mehrere Integratoren mit; in der Praxis wird fast immer CVODE aus dem SUNDIALS-Paket benutzt — ein implizites, adaptives Verfahren auf Basis von BDF (Backward Differentiation Formulas) der Ordnung 1 bis 5. Explizite Verfahren (Explicit Euler, Runge-Kutta 45) sind als Referenz- und Debugging-Integratoren ebenfalls vorhanden.

CVODE schätzt nach jedem Versuch den lokalen Fehler ab und vergleicht ihn mit einer vorgegebenen relativen und absoluten Fehlertoleranz. Ist der Fehler zu groß, wird der Schritt verworfen und mit kleinerem dt (oder niedrigerer Ordnung) wiederholt; ist er bequem klein, dürfen Schrittweite und Ordnung wachsen. So findet der Solver automatisch den Kompromiss aus Genauigkeit und Geschwindigkeit.

Explizit vs. implizit — und warum Gebäude implizit gerechnet werden

Ob ein Zeitschritt eine Newton-Iteration und ein lineares Gleichungssystem braucht, hängt allein davon ab, wo das unbekannte y_{n+1} in der Formel auftaucht.

Explizit — nur alte, bekannte Werte rechts. Man setzt ein und rechnet aus, eine einzige ydot-Auswertung pro Schritt, nichts aufzulösen:

$y_{n+1} = y_n + dt \cdot \texttt{ydot}(t_n,\, y_n)$

Implizit — das Unbekannte steht auch rechts und muss aufgelöst werden:

$y_{n+1} = y_n + dt \cdot \texttt{ydot}(t_{n+1},\, y_{n+1})$

Wenn explizit pro Schritt billiger ist — warum dann implizit? Wegen Steifigkeit. Gebäudemodelle vereinen sehr schnelle (Luft, dünne Schichten) und sehr langsame (massive Wände, Estrich) Zeitkonstanten. Explizite Verfahren müssen dann aus Stabilitätsgründen winzige Schritte machen — die Schrittweite wird nicht von der gewünschten Genauigkeit bestimmt, sondern von der schnellsten Zeitkonstante. Tausende billige Mini-Schritte sind am Ende deutlich teurer als die wenigen großen, genauigkeitsgetriebenen Schritte eines impliziten Verfahrens.

Den Effekt sehen Sie unmittelbar in den ersten beiden Demos: Schieben Sie dt über die Stabilitätsgrenze, beginnt die explizite Lösung zu oszillieren und explodiert — die implizite bleibt ruhig. In der Wand-Demo macht ein feineres Netz (großes N) das Problem zusätzlich steifer, und die explizite Fourier-Zahl $r = \alpha\,dt/\Delta x^2$ sprengt schnell die Grenze $r \le \tfrac12$.

Drei geschachtelte Schleifen

Das ist die wichtigste mentale Landkarte. Der implizite Solver besteht aus drei Ebenen, die ineinanderstecken — genau die drei Schleifen aus dem interaktiven Schema oben:

• Außen — Zeitschleife (CVODE): wählt die Schrittweite dt und geht von t nach t+dt.
• Mitte — Newton-Iteration: löst die nichtlineare implizite Gleichung nach y_{n+1} (meist 1–3 Iterationen).
• Innen — Lineares Gleichungssystem (LES): in jeder Newton-Iteration wird einmal $J \cdot \Delta y = -F$ gelöst.

Die Newton-Iteration

Die implizite Gleichung schreibt man als „= 0”:

$F(y) = y - y_n - dt \cdot \texttt{ydot}(y) = 0$

Newton verbessert eine Schätzung schrittweise, indem es lokal linearisiert:

$J \cdot \Delta y = -F(y), \qquad \text{dann:}\quad y \leftarrow y + \Delta y$

Dabei ist $J$ die Jacobi-Matrix — die Ableitung von $F$ nach $y$, also „wie reagiert jede Bilanz auf eine kleine Änderung jeder Zustandsgröße”. Diese Matrixgleichung ist genau das lineare Gleichungssystem — und $J$ ist dessen Koeffizientenmatrix (das $A$ in $A\,\Delta y = b$), nicht etwa ein Hilfsmittel zur Fehlerschätzung. Die adaptive Schrittweitensteuerung schätzt den Zeitschritt-Fehler getrennt davon über die Toleranzen ab und benötigt $J$ dafür nicht. Bei $n$ Zuständen ist $J$ eine $n\times n$-Matrix, $-F$ und $\Delta y$ sind Vektoren der Länge $n$ — ein quadratisches System aus $n$ Gleichungen.

Wo ydot einfließt. Beide Bestandteile dieses Systems werden aus ydot gebaut — ydot ist kein Nebenstrang. Die rechte Seite ist direkt das negative Residuum, $-F = -\bigl(y - y_n - dt\,\dot y\bigr)$, also eine ydot-Auswertung am aktuellen Rateversuch. Die Matrix ist

$J = \frac{\partial F}{\partial y} = I - dt\,\frac{\partial \dot y}{\partial y}$

und die Empfindlichkeit $\partial \dot y/\partial y$ entsteht ihrerseits aus ydot-Auswertungen: NANDRAD stört jede Zustandsgröße ein wenig, lässt die Head→Mitte→Tail-Pipeline mit dem gestörten Zustand erneut laufen und bildet den Differenzenquotienten — das ergibt eine Spalte von $J$. Jede Newton-Iteration ruft ydot also mehrfach auf: einmal für die rechte Seite und zusätzlich für die Spalten der Jacobi-Matrix (beschleunigt durch Sparsity und Graph-Coloring, siehe unten).

Das lineare Gleichungssystem (LES)

In jeder Newton-Iteration muss $J \cdot \Delta y = -F$ gelöst werden. $J$ ist eine $n \times n$-Matrix; bei tausenden Unbekannten ist das der rechenintensivste Teil der ganzen Simulation.

Woher kommt die Jacobi-Matrix?

Meist wird $J$ numerisch über finite Differenzen gebildet: Man stört jede Zustandsgröße ein wenig, schaut, wie sich $F$ ändert, und erhält so spaltenweise die Matrix. Bei tausenden Spalten wäre das naiv sehr teuer — der Trick ist Sparsity (Dünnbesetztheit): Eine Wandschicht koppelt nur mit ihren direkten Nachbarn, ein Raum nur mit seinen angrenzenden Bauteilen. Die allermeisten Einträge von $J$ sind also Null. NANDRAD kennt das Besetzungsmuster vorab und stört per Graph-Coloring viele unabhängige Spalten gleichzeitig — statt tausender braucht es oft nur eine Handvoll Auswertungen.

Die tridiagonale Struktur einer eindimensionalen Wand-Diskretisierung ist der einfachste Fall dieser Dünnbesetztheit. In der Wand-Demo wird genau dieses tridiagonale System gezeigt und Zelle für Zelle aufgeschlüsselt — die dort dargestellte Matrix mit $1+2r$ auf der Diagonale und $-r$ daneben ist die Jacobi-Matrix, also die Koeffizientenmatrix des LES.

Direkte vs. iterative Löser

Es gibt zwei grundlegend verschiedene Arten, $J \cdot \Delta y = -F$ zu lösen:

• Direkte Löser berechnen die exakte Lösung über eine Faktorisierung (LU) — gut für kleine, volle ebenso wie für große, dünnbesetzte Systeme.
• Iterative Löser nähern die Lösung schrittweise an (Krylov-Verfahren); sie brauchen einen Präkonditionierer, der das System vorab „glättet”.

Lässt man diese Felder im Projekt leer, wählt NANDRAD sinnvolle Voreinstellungen (in der Regel einen direkten Sparse-Löser, bei großen Netzwerken einen iterativen Löser mit Präkonditionierer). Für die meisten Anwendungen ist das die richtige Wahl — an diesen Parametern dreht man nur, wenn eine Simulation zu langsam ist oder nicht konvergiert.

Vom Zustand zur Lösung in fünf Sätzen

1. Der Solver bilanziert Energie pro Raum und pro Wandschicht; daraus ergibt sich ein großes System gewöhnlicher Differentialgleichungen $dy/dt = \texttt{ydot}(t, y)$.
2. Die Funktion ydot wird in drei Stufen ausgewertet (Energie → Temperatur, Flüsse, Summe → ydot); die Reihenfolge bestimmt ein automatischer Abhängigkeitsgraph.
3. Integriert wird mit CVODE (implizites BDF-Verfahren) und adaptiver Schrittweite plus Fehlerkontrolle — nötig, weil Gebäudemodelle steif sind.
4. Jeder implizite Zeitschritt erfordert eine Newton-Iteration, und jede Newton-Iteration ein lineares Gleichungssystem $J \cdot \Delta y = -F$ (die drei geschachtelten Schleifen).
5. Das LES wird je nach Größe und Dünnbesetztheit direkt (Dense/KLU) oder iterativ (GMRES/BiCGStab + ILU) gelöst; die Jacobi-Matrix entsteht numerisch und nutzt die Sparsity per Graph-Coloring aus.

Was das für die Praxis bedeutet

Für die tägliche Arbeit im Ingenieurbüro müssen Sie diese Maschinerie nicht im Detail bedienen — aber es hilft zu verstehen, warum eine validierte dynamische Simulation belastbare Ergebnisse liefert und an welchen Stellschrauben man im Zweifel dreht. Der NandradSolver ist der Rechenkern, auf dem VICUS Buildings aufsetzt; die zugrundeliegende Methodik ist quelloffen über das Forschungsprojekt SIM-VICUS dokumentiert und in über 20 Jahren an der TU Dresden entwickelt worden. Wer tiefer einsteigen will, findet weiterführende Grundlagen in unseren Artikeln zur dynamischen Gebäudesimulation, zum Unterschied stationär vs. dynamisch und zu den Validierungsstandards.