Simulation and optimisation of district heating and cooling networks requires efficient and realistic models of the individual network elements in order to correctly represent heat losses or gains, temperature propagation and pressure drops. Due to more recent thermal networks incorporating meshing decentralised heat and cold sources, the system often has to deal with variable temperatures and mass flow rates, with flow reversal occurring more frequently. This paper presents the mathematical derivation and software implementation in Modelica of a thermo-hydraulic model for thermal networks that meets the above requirements and compares it to both experimental data and a commonly used model. Good correspondence between experimental data from a controlled test set-up and simulations using the presented model was found. Compared to measurement data from a real district heating network, the simulation results led to a larger error than in the controlled test set-up, but the general trend is still approximated closely and the model yields results similar to a pipe model from the Modelica Standard Library. However, the presented model simulates 1.7 (for low number of volumes) to 68 (for highly discretized pipes) times faster than a conventional model for a realistic test case. A working implementation of the presented model is made openly available within the IBPSA Modelica Library. The model is robust in the sense that grid size and time step do not need to be adapted to the flow rate, as is the case in finite volume models.