Every thermo-hydraulic simulation stands or falls with its computational core, the piece of software that turns network topology, pipe data, load profiles and control logic into a physically reliable pressure, mass-flow and temperature history. In VICUS Districts this core is the VICUS Solver, a simulation solver matured over more than 20 years of research at TU Dresden — the same core that powers VICUS Buildings. This article lifts the hood for the district heating case: step by step it shows how the solver describes a network as a system of equations and integrates it through time.

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

The core idea in one sentence

The solver describes the heat network as a graph of nodes and edges and solves two tightly coupled subproblems on it: the hydraulics — which pressures and mass flows establish themselves? — and the heat transport — how do the temperatures travel through the network? The heart of it is always the same question:

„How fast is the stored energy changing in every fluid volume of the network?”

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 network as a graph

The basis of every calculation is the network topology as a graph: nodes are branching and connection points, edges are the network elements — pipes, pumps, valves, the heat exchangers of the transfer stations, producers. Every element brings two models: a hydraulic one (which pressure difference establishes itself at which mass flow?) and a thermal one (what happens to the temperature of the fluid passing through?).

This separation is not a side note but the central architectural decision of the solver — because the two subproblems live on fundamentally different time scales.

Two subproblems: fast hydraulics, sluggish thermal transport

The hydraulics: steady-state and algebraic

Pressure waves propagate through nearly incompressible water at the speed of sound. From the perspective of the slowly travelling temperatures, the flow field therefore establishes itself instantly. That is why the hydraulics are not integrated over time but solved as a steady-state, algebraic system of equations in every evaluation:

• Mass conservation at every node $i$ (Kirchhoff’s first law): $\sum_k \dot{m}_{i,k} = 0$
• Pressure loss per element, e.g. for pipes according to Colebrook-White — depending on mass flow and temperature, since the density and viscosity of the fluid change with it (see pressure loss calculation)

This nonlinear system is solved with the Newton-Raphson method. Bidirectional flows — the flow direction reverses, for instance with decentralised feed-in — are automatically included, because the pressure-loss equation is formulated sign-correctly.

The thermal transport: dynamic and time-dependent

Temperatures change comparatively slowly: the fluid stores energy, and a temperature wave takes minutes to hours to travel through the network, depending on flow velocity. Heat transport is therefore formulated as a system of ordinary differential equations — the energy balance of every fluid volume $k$:

$C_k \cdot \frac{dT_k}{dt} = \dot{H}_{\text{in},k} - \dot{H}_{\text{out},k} - \dot{Q}_{\text{loss},k}$

with the enthalpy flows $\dot H = \dot m \, c \, T$ and the heat loss to the ground $\dot{Q}_{\text{loss},k} = UA \cdot (T_k - T_\text{ground})$.

Why not put everything into one single large system of equations? That is the route of generic solvers (e.g. in Modelica): one monolithic Jacobian with many direct dependencies. The VICUS Solver takes the considerably more performant second route: the hydraulic system is solved as a separate subproblem within the thermal calculation — two smaller, specialised systems of equations instead of one large one.

Energy instead of temperature: the conserved quantity

As in the building case, the solver does not balance temperature directly but the stored energy — because energy is a conserved quantity whose inflows and outflows can be summed unambiguously. Temperature is a quantity derived from it.

What is in the state vector y?

The dynamic pipe model splits every pipe into segments so that the propagation of temperature waves is captured correctly; each segment is a separate entry in y. The thermal capacity of the pipe wall is attributed to the fluid via a lumped model — the inertia is correct without an additional computational grid for the pipe wall. For large networks this quickly amounts to thousands of unknowns.

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

The recipe is the same as for every conservation law:

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 pipe segment

A segment $k$ with fluid mass $M$ and heat capacity $c$ is passed by the mass flow $\dot m$ and loses heat to the ground:

$M\,c\,\frac{dT_k}{dt} = \dot m\,c\,\bigl(T_{k-1} - T_k\bigr) \;-\; UA\,\bigl(T_k - T_\text{ground}\bigr)$

On the left the change in storage, on the right advection (the warm fluid from the previous segment) minus loss. Solved for $dT_k/dt$, this is the ydot for this segment — derived from energy conservation plus flux laws, not invented. Chaining many segments together produces exactly the model that runs in the live demo below: a temperature front travelling along a pipe split into 50 segments.

The runtime evaluation — hydraulics included

There is no closed ydot formula sheet in the network case either. The integrator hands over a y and asks for the rate of change; internally the familiar three-stage pipeline runs:

The middle stage is remarkable: every single ydot evaluation contains a complete hydraulic network calculation, because without current mass flows there are no enthalpy flows. A dependency graph automatically sorts the sub-models into the correct order — controllers that react to quantities they influence themselves are grouped together and solved iteratively.

The schema, interactively

The visualisation below shows the functional schema of the VICUS Solver with its three nested loops — in the network case it is the same as in the building case, because it is the same computational core. The live demo underneath shows the building block of every network model: a cold pipe, split into segments, into which hot supply water flows from $t = 0$. You can watch the temperature front travel along the pipe at the flow velocity while heat losses damp it on the way — and see where the stability limit of explicit methods lies.

Time integration: CVODE and the stiffness of the network

Integration uses CVODE from the SUNDIALS package: an implicit, adaptive multi-step method based on BDF (Backward Differentiation Formulas) of order 1 to 5. The step size is not constant — during calm summer operation the solver takes large steps, during the morning load ramp-up or when a producer comes online it takes small ones. The output (hourly values, say) is decoupled from this and cleanly interpolated.

Why implicit? For the same reason as with buildings: stiffness. A short, fast-flowing pipe segment has a time constant of seconds; the ground around the route responds over weeks. Explicit methods would have to follow the fastest time constant and take tiny steps — an annual simulation would become unusably slow. The implicit method stays stable with large steps; the step size is governed solely by the desired accuracy.

Three nested loops — plus a subproblem

The mental map of the implicit solver applies unchanged:

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

The network case adds the special feature that every ydot evaluation in turn solves the hydraulic subproblem via Newton-Raphson. There are thus two separate Jacobians, each with few dependencies — the hydraulic one is formed analytically for simple elements and numerically for controlled components.

Both systems are sparse: a pipe segment couples thermally only with its neighbours along the flow path, a node hydraulically only with the connected edges. The solver exploits this via sparse solvers; for very large networks, iterative Krylov methods with a preconditioner are used (GMRES + ILU). For most projects the solver picks these settings itself — you only touch them when a simulation is too slow or fails to converge.

From network model to solution in five sentences

1. The solver describes the heat network as a graph; every element brings a hydraulic and a thermal model.
2. The hydraulics (pressures, mass flows) are a steady-state algebraic system solved by Newton-Raphson in every evaluation; the thermal transport (temperatures) is a system of ordinary differential equations $dy/dt = \texttt{ydot}(t, y)$.
3. ydot balances the energy of every fluid volume — advection minus heat loss — and is evaluated at runtime in three stages, including the hydraulic subproblem.
4. Integration uses CVODE (implicit BDF method) with adaptive step size plus error control — necessary because network models are stiff.
5. Each implicit time step requires a Newton iteration, each Newton iteration a sparse linear system of equations, solved directly (sparse LU) or iteratively (GMRES + ILU).

What this means in practice

For day-to-day network planning you do not need to operate this machinery in detail — but it explains why a genuine thermo-hydraulic annual simulation needs computing time and is nevertheless finished in minutes rather than hours: the domain-specific split into hydraulic and thermal subproblems is the performance lever that makes even large networks with hundreds of participants, decentralised feed-in and prosumer structures reliably computable. The VICUS Solver is the computational core that VICUS Districts is built on; the underlying methodology is documented as open source through the SIM-VICUS research project.

To dig deeper: the article on thermo-hydraulic simulation places the methodology within network planning, pressure loss calculation shows the hydraulic equations in detail, and VICUS Solver: the Engine Behind VICUS Buildings explains the same core from the building perspective — with further interactive demos on the Newton iteration and the linear system of equations.