The Wanda software package developed by Deltares can used for simulating both steady state and transient fluid flow in pipeline systems. Steady state simulations are used for initial system design and transient flow simulations are used for doing water hammer analysis in pipeline systems. For both types of simulations a system consisting of both linear and non-linear equations needs to be solved for the main unknown quantities, flowrate and head. This system is solved by linearising the equations using theNewton-Raphson method and solving the resulting system of linear equations. Currently, this is done by using a matrix solver from the proprietary IMSL numerical library which requires a paid license. The problem is that this solver sometimes either crashes or gets stuck in an infinite loop when dealing with singular matrices, while the proprietary nature of the library only allows for limited troubleshooting. The solution method therefore requires a redesign which should improve its robustness and maintainability. On the other hand, no ground should be yielded in terms of solution accuracy and performance.

The singularity of the matrices are caused by quantities being underdetermined either due to user error in network design or phase changes such as a closing valve. In this report, a graph-theoretic approach is taken to detect these structural singularities in the form of determining a maximum size matching in a graph representing the system of equations. This approach gives information about which quantity is undetermined where in the pipeline system. As an alternative, condition number estimation is implemented. Furthermore, the IMSL library is replaced by LAPACK, which is a lightweight and versatile alternative. The open source nature of LAPACK and its permissive license ensure its maintainability. Since thematrices are banded, the band version of the LAPACK algorithms can be used. The new solutionmethod is compared to IMSL and evaluated in terms of robustness, accuracy and performance.

The graph-theoretic method resolves the robustness issues of the IMSL-based method, while showing great performance. The information it gives is used to either correct the matrix and continue the simulation or output an appropriate error message, ensuring a user-friendly experience. Condition number estimation is too slow while also not being useful for further matrix correction purposes and is therefore disregarded. To improve the accuracy of LAPACK, iterative refinement is used. The maximum relative error measured over all the test cases was about 2.5%, resulting from an ill-conditioned case. Overall, the accuracy compared to IMSL is good, so that users will not be able to notice large difference in solution quality between the two solution methods. To improve the performance of LAPACK, Reverse Cuthill-Mckee (RCM) is applied to reduce the bandwidth of the matrices. Using several test cases, the LAPACK performance is shown to be similar to that of IMSL when using RCM. In that respect, the transition from IMSL to LAPACK should also be flawless. A vendor-optimised LAPACK implementation did not yield any significant performance gains. It can be concluded that the new solution method is an improvement in terms of robustness and maintainability, while showing similar solution accuracy and performance.

Future work can focus on improving the performance by calculating the RCM permutation only once every time step, as well as keeping the Jacobian constant during each time step. As matrix bandwidth determines the performance of the matrix solver, better methods to reduce the matrix bandwidth could also yield significant improvements. This could be achieved by either ordering the components in the pipeline system a priori, or simply using matrix bandwidth reduction techniques. Optimised LAPACK libraries as well as a different numerical library such asMUMPS could improve the performance even further.