Every dynamic building simulation stands or falls with its computational core — the piece of software that turns geometry, building components, weather and use into a physically reliable temperature and load profile. In VICUS Buildings this core is NANDRAD, a simulation solver matured over more than 20 years of research at TU Dresden. This article lifts the hood: step by step it shows what the NandradSolver computes and how it does it — from the physical balance equation via time integration to the linear system of equations at its heart.

You do not need a numerics degree for this. The central idea fits in a single sentence and can be „cranked through” yourself with the interactive demos below.

The core idea in one sentence

The solver describes the building (and possibly a pipe network) as one large system of coupled ordinary differential equations and integrates this system step by step through time. The heart of it is always the same question:

„How fast is the stored energy changing in every room and in every wall layer?”

Once you know this rate of change, you can take a small step into the future, compute the new temperatures and start over. In vector notation — exactly how the time integrator „sees” it — this reads:

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

Here $\mathbf{y}$ is the state vector (a long list of all conserved quantities), $\texttt{ydot}$ their rate of change and $t$ the time.

The schema and three live demos

The interactive visualisation below summarises the whole process — and makes it tangible. At the top the functional schema with the three nested loops (hover the boxes to see the relevant source files; toggle between implicit and explicit). Beneath it three live demos in which you can turn the time step, grid fineness and heat-exchanger parameters yourself and find the stability limit.

The rest of the article explains the individual building blocks in detail. You can scroll back up at any time and try out what you have learned.

Energy instead of temperature: the conserved quantity

A central design decision in NANDRAD: the solver balances energy (more precisely the stored internal energy), not temperature directly. The reason is physically clean — energy is a conserved quantity. All energy fluxes flowing into and out of a room can be summed unambiguously. Temperature is then a derived quantity, recovered from the energy and the heat capacity.

The general form of every balance equation is:

$\frac{d(\text{stored energy})}{dt} = \sum \text{inflows} - \sum \text{outflows}$

What is in the state vector y?

y is simply everything concatenated that changes over time and stores energy:

A single wall is therefore split internally into many thin layers — each layer is a separate entry in y. For a larger building this quickly amounts to several thousand unknowns. You can generate this „one wall into many coupled equations” decomposition yourself in the second demo (heat conduction through a wall) with the N slider.

Where does ydot come from? The balance, not a guess

A fair question: $dy/dt = \texttt{ydot}(t, y)$ is stated — but where does the right-hand side come from? The answer: ydot is not guessed, it is derived from a conservation law. The recipe is always the same:

1. Choose the conserved quantity y (energy, mass, …) — not the temperature.
2. Write the balance: d(stored)/dt = Σ inflows − Σ outflows.
3. Express each flux through the state (physical laws).
4. Solve for dy/dt — that is ydot.

Example: a room cools down

The stored quantity is the internal energy $U = C\,T$ with heat capacity $C$. It only changes through the heat flow $Q$ across the envelope, which is proportional to the temperature difference:

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

With $\frac{dU}{dt} = C\,\frac{dT}{dt}$, solving for the rate gives:

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

So ydot is derived, not invented — energy conservation plus one flux law. Exactly this single equation runs in the first demo („a room cools down”): you see the integrator work its way along the tangent through time, compared against the analytic solution.

And for a real building?

The crucial insight: in NANDRAD there is no single closed ydot formula sheet. With thousands of states a hand-written formula would be impossible. Instead the balance is the same for every state, but the sum of fluxes is assembled at runtime — every sub-model throws its contribution into the pot: conduction between layers, solar gains through windows, ventilation and infiltration, internal gains, heating and cooling power, pipe flow. ydot is therefore a function you do not write down but evaluate.

The three evaluation stages: Head → Middle → Tail

The integrator hands the model a y and asks for the corresponding ydot. Internally this happens in three stages — exactly the balance recipe, cast into code:

In the schema above this data flow is shown as a three-stage pipeline inside the Newton loop.

Order of the models: the dependency graph

There are many small sub-models, and they depend on one another: the window model needs the outdoor temperature from the climate model, the room balance needs the heating power from the HVAC model. In what order must one compute?

NANDRAD solves this with a dependency graph. Each model declares which quantities it needs as input and which it produces as a result. From this a correct evaluation order is computed automatically (topological sort). Two special cases are noteworthy:

• Parallelisation: models that are independent of one another sit in the same „level” of the graph and are evaluated in parallel (e.g. many walls at once).
• Cycles: if models depend on each other mutually (such as a controller reacting to a quantity it influences itself), they are grouped together and solved iteratively as a group.

Time integration: taking a step into the future

Now we have a function that yields ydot for any y. The job of time integration is to reconstruct the time history from it — the path from the slope. In an outer loop the solver advances step by step from t to t+dt, commits each accepted step and writes out any due results.

Importantly, the step size dt is not constant. The integrator chooses it itself — large steps at night when little changes, very small ones when the heating kicks in. The output (hourly values, say) is decoupled from this: between the actual computation steps, values are cleanly interpolated onto the requested output times.

Which integrator? CVODE as the default

NANDRAD ships several integrators; in practice CVODE from the SUNDIALS package is used almost always — an implicit, adaptive method based on BDF (Backward Differentiation Formulas) of order 1 to 5. Explicit methods (Explicit Euler, Runge-Kutta 45) are also available as reference and debugging integrators.

After each attempt CVODE estimates the local error and compares it with a prescribed relative and absolute error tolerance. If the error is too large, the step is rejected and retried with a smaller dt (or lower order); if it is comfortably small, step size and order may grow. This is how the solver automatically finds the compromise between accuracy and speed.

Explicit vs. implicit — and why buildings are computed implicitly

Whether a time step needs a Newton iteration and a linear system of equations depends solely on where the unknown y_{n+1} appears in the formula.

Explicit — only old, known values on the right. You substitute and compute, a single ydot evaluation per step, nothing to solve:

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

Implicit — the unknown also appears on the right and must be solved for:

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

If explicit is cheaper per step — why implicit? Because of stiffness. Building models combine very fast (air, thin layers) and very slow (massive walls, screed) time constants. Explicit methods then have to take tiny steps for stability reasons — the step size is governed not by the desired accuracy but by the fastest time constant. Thousands of cheap mini-steps end up far more expensive than the few large, accuracy-driven steps of an implicit method.

You can see the effect immediately in the first two demos: push dt past the stability limit and the explicit solution starts to oscillate and explodes — the implicit one stays calm. In the wall demo a finer grid (large N) makes the problem even stiffer, and the explicit Fourier number $r = \alpha\,dt/\Delta x^2$ quickly bursts the limit $r \le \tfrac12$.

Three nested loops

This is the most important mental map. The implicit solver consists of three levels nested inside each other — exactly the three loops from the interactive schema above:

• Outer — time loop (CVODE): picks the step size dt and goes from t to t+dt.
• Middle — Newton iteration: solves the nonlinear implicit equation for y_{n+1} (usually 1–3 iterations).
• Inner — linear system of equations (LSE): in every Newton iteration, $J \cdot \Delta y = -F$ is solved once.

The Newton iteration

The implicit equation is written as „= 0”:

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

Newton improves a guess step by step by linearising locally:

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

Here $J$ is the Jacobian — the derivative of $F$ with respect to $y$, i.e. „how does each balance react to a small change of each state quantity”. This matrix equation is exactly the linear system of equations — and $J$ is its coefficient matrix (the $A$ in $A\,\Delta y = b$), not a tool for error estimation. Adaptive step-size control estimates the time-step error separately via the tolerances and does not need $J$ for that. With $n$ states, $J$ is an $n\times n$ matrix while $-F$ and $\Delta y$ are vectors of length $n$ — a square system of $n$ equations.

Where ydot enters. Both parts of this system are built from ydot — it is not a side branch. The right-hand side is directly the negative residual, $-F = -\bigl(y - y_n - dt\,\dot y\bigr)$, i.e. one ydot evaluation at the current guess. The matrix is

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

and the sensitivity $\partial \dot y/\partial y$ is itself obtained from ydot evaluations: NANDRAD perturbs each state quantity slightly, re-runs the Head→Middle→Tail pipeline with the perturbed state and forms the finite difference — that yields one column of $J$. So each Newton iteration calls ydot several times: once for the right-hand side and additionally for the columns of the Jacobian (accelerated by sparsity and graph colouring, see below).

The linear system of equations (LSE)

In every Newton iteration $J \cdot \Delta y = -F$ must be solved. $J$ is an $n \times n$ matrix; with thousands of unknowns this is the most compute-intensive part of the whole simulation.

Where does the Jacobian come from?

Usually $J$ is formed numerically via finite differences: each state quantity is perturbed slightly, the change in $F$ is observed, and the matrix is obtained column by column. With thousands of columns this would naively be very expensive — the trick is sparsity: a wall layer couples only with its direct neighbours, a room only with its adjacent components. The vast majority of entries in $J$ are therefore zero. NANDRAD knows the sparsity pattern in advance and, via graph colouring, perturbs many independent columns simultaneously — instead of thousands it often needs only a handful of evaluations.

The tridiagonal structure of a one-dimensional wall discretisation is the simplest case of this sparsity. The wall demo shows exactly this tridiagonal system and breaks it down cell by cell — the matrix shown there, with $1+2r$ on the diagonal and $-r$ next to it, is the Jacobian, i.e. the coefficient matrix of the LSE.

Direct vs. iterative solvers

There are two fundamentally different ways to solve $J \cdot \Delta y = -F$:

• Direct solvers compute the exact solution via a factorisation (LU) — good for small, full as well as large, sparse systems.
• Iterative solvers approximate the solution step by step (Krylov methods); they need a preconditioner that „smooths” the system beforehand.

If you leave these fields empty in the project, NANDRAD picks sensible defaults (usually a direct sparse solver, or an iterative solver with a preconditioner for large networks). For most applications this is the right choice — you only touch these parameters when a simulation is too slow or fails to converge.

From state to solution in five sentences

1. The solver balances energy per room and per wall layer; this yields a large system of ordinary differential equations $dy/dt = \texttt{ydot}(t, y)$.
2. The function ydot is evaluated in three stages (energy → temperature, fluxes, sum → ydot); the order is determined by an automatic dependency graph.
3. Integration uses CVODE (implicit BDF method) with adaptive step size plus error control — necessary because building models are stiff.
4. Each implicit time step requires a Newton iteration, and each Newton iteration a linear system of equations $J \cdot \Delta y = -F$ (the three nested loops).
5. Depending on size and sparsity, the LSE is solved directly (Dense/KLU) or iteratively (GMRES/BiCGStab + ILU); the Jacobian is formed numerically and exploits the sparsity via graph colouring.

What this means in practice

For day-to-day work in an engineering office you do not need to operate this machinery in detail — but it helps to understand why a validated dynamic simulation delivers reliable results and which screws to turn in case of doubt. The NandradSolver is the computational core that VICUS Buildings builds on; the underlying methodology is documented as open source through the SIM-VICUS research project and was developed over more than 20 years at TU Dresden. To dig deeper, see our articles on dynamic building simulation, on the difference between steady-state and dynamic, and on validation standards.