## Abstract

This paper presents the development of an efficient, interactive and fully extendable computer code to analyse hydraulic transients – water hammer – using a generic object-oriented modelling and simulation (OOMS) tool to simplify and shorten the development process. OOMS tools provide continuous simulation features by means of ordinary differential equation (ODE) and differential algebraic equation (DAE) solvers, and discrete simulation features by means of event handling algorithms. Previous OOMS tool applications to simulate fluid flow in pipe networks primarily used semi-discrete methods with continuous solvers. A novel aspect of this work is the application of the method of characteristics (MOC) using the discrete simulation features of the OOMS tool. When compared to a code developed with an OOMS tool using a semi-discrete method and to commercial software that uses the MOC, the new code shows much higher accuracy and performance than the former and is similar to the latter in accuracy and calculation time. It has a graphical user interface, and its modularity makes it easily extendable with new components and algorithms. An academic version is available on github.

## HIGHLIGHTS

Physical modeling languages (PMLs), such as Modelica and EcosimPro, shorten the development of simulation software.

PMLs rely on generic solvers, which are not efficient for specific problems (e.g. hydraulic transients).

The purpose of the research is to enable the implementation of

*ad hoc*solvers using PMLs.Results show the feasibility of the new approach and its applicability to hydraulic transients

### Graphical Abstract

## NOMENCLATURE

pipe cross-sectional area (m

^{2})valve exit area when fully open (m

^{2})wave speed of the water hammer waves (m/s)

pipe compressibility constant (s/m

^{2})quantity for negative characteristics used for 1D MOC solution (s/m

^{2})quantity for positive characteristics used for 1D MOC solution (s/m

^{2})quantity for negative characteristics used for 1D MOC solution (m)

quantity for positive characteristics used for 1D MOC solution (m)

discharge coefficient in exit valve (unitless)

Courant number (unitless)

pipe internal diameter (m)

Darcy friction factor (unitless)

gravitational acceleration (m/s

^{2})piezometric head (m)

vapour pressure head (m)

friction decay coefficient for the calculation of the transient friction (unitless)

user-defined constant for the calculation of the artificial viscosity term (unitless)

pipe length (m)

static pressure (N/m

^{2})volume flow (m

^{3}/s)upstream volume flow (m

^{3}/s)downstream volume flow (m

^{3}/s)resistance coefficient (s

^{2}/m^{5})distance along the pipe (m)

time (s)

vapour cavity volume (m

^{3})artificial viscosity head (m)

elevation (m)

linearization constant for the friction term (unitless)

loss coefficient (unitless)

permissible variation in the wave speed of pipe

*J*(unitless)density of liquid (kg/m

^{3})relative valve opening (unitless)

**Subscripts**

## INTRODUCTION

Transient fluid flow in pipes can generate pressure waves that travel through the body of the fluid at the local speed of sound. Pressure surges and other related phenomena are referred to as water hammer when the fluid is a liquid. Water hammer can occur in a wide variety of situations, such as in aerospace propulsion systems, hydro-power stations, power plant steam generation/cooling systems, district heating systems, industrial plants, water supply systems, irrigation networks, domestic plumbing, and oil and chemical transport pipelines.

Nowadays, there are many commercial software tools available on the market for calculating water hammer. These tools have been developed using general-purpose programming languages (GPPLs) such as Fortran, C, C ++ and Java. GPPLs are the tools of choice for scientific software for specific fields because they provide the flexibility and the ability to implement specific domain solvers that are efficient, accurate and robust. Their main drawback is that they require enormous amounts of time to write, document and debug the code. This problem has been exacerbated by the emergence of graphical user interfaces (GUIs). Although GUIs make the software easier to use and learn, they have significantly increased the total effort required to create scientific applications.

Object-oriented modelling and simulation (OOMS) tools such as EcosimPro (EAI 2015) and the Modelica-based tools (Mattsson & Elmqvist 1997) were created with the goal of reducing the effort required to develop simulation applications. These tools are designed around equation-based languages, conceived to represent models by differential algebraic equations (DAEs), difference equations and discrete events.

The traditional approach to simulation codes using OOMS tools is to formulate the problem as a DAE system. This study uses a new approach: instead of formulating the simulation problem using DAEs, the discrete simulation features of the OOMS tool were used to implement solvers specifically designed for the given simulation domain. This mixed approach combines the advantages of OOMS tools, i.e., shortened and simplified development, and the advantages of GPPLs, i.e., to provide the tool with specific solvers, which are efficient and accurate for the given field. The development of a water hammer code, named LiqHammer, using the OOMS tool EcosimPro demonstrates the feasibility of the approach. This paper outlines the design of the LiqHammer code, which implements the MOC and the Algebraic Water Hammer (AWH) method (Wylie *et al.* 1993), and presents the results obtained for a few example cases. A comparison of the new code with other water hammer codes shows that it is as efficient as the codes developed using GPPLs, and the development effort is significantly reduced due to the use of an OOMS tool.

The rest of this article is structured as follows. Next section provides a brief overview of OOMS. Then, after briefly describing the state of the art of water hammer software, the MOC for water hammer analysis is presented, there is nothing new in this presentation and has been included only for the sake of completeness. Following section presents an overview of the new implementation of the MOC in an OOMS tool from the points of view of both developers and users. Afterwards, the comparison of the results of the new approach with the values of the bibliography references and the results calculated by commercial water hammer analysis codes are shown and discussed. The last section concludes the article by summarizing the findings of the study.

## OBJECT-ORIENTED MODELLING AND SIMULATION

OOMS tools typically offer a modelling language, which is object-oriented, and a simulation environment with two clearly differentiated parts: solvers and algorithms for finding good solutions to the mathematical simulation problem, and a graphical user interface with interactive tools for building models, running experiments and analysing results.

### Modelling language

The main idea behind the OOMS tools is to support the creation of libraries of reusable unit models, usually named component types or classes, by means of declarative modelling languages. Complex models can then be built by instantiating and connecting components in a process known as aggregation.

OOMS tools can be tailored to simulate specific domains. Components usually represent physical parts of the modelled system. For example, a library for hydraulic simulation includes components such as pipes, valves, junctions, tanks and pumps.

Components have ‘connection ports’ that are endpoints for different types of connections, such as fluid flow connections, heat transfer connections, electrical current connections (wires), mechanical links and control signals. A system model is built at the topological level (or flowsheet level) by connecting components in the same way as they are physically connected.

### Mathematical formulation of OOMS tools

#### Differential algebraic equations

*t*is time. DAEs consist of differential equations, (the vector) and algebraic constraints (the vector) and are a generalization of ordinary differential equations (ODEs). They are used to represent the continuous dynamic behaviour of the system. Numerical methods for integration of DAEs appeared in the 1970s, and many different DAE solvers have been developed over the last decades.

#### Events and discontinuous behaviour

Events and their associated actions define the discrete behaviour of the system. Events have zero duration and consist of a trigger condition and the action to be performed when the action is triggered. The trigger condition is a Boolean expression and the event occurs when the condition changes from false to true. Depending on the trigger condition, events can be classified into *time-events* and *state-events*.

A *time-event* is one whose trigger condition only depends on the time variable; therefore, its time of occurrence is known in advance, and the end time of the integration step can be set equal to the time of the event. In contrast, the trigger condition of a *state-event* depends on any continuous-time variable (), so their time of occurrence cannot be predicted in advance. The trigger conditions of *state-events* are monitored at the end of the integration steps of the DAE solver. If any condition changes to true, the tool finds the precise time of the event and adjusts the end time of the last integration step to it.

The event action describes how the new and old state relate to each other. Events in OOMS tools are specified using *when-clauses*, *conditional equations* and *delayed assignments*.

*when-clauses,*the value of the discrete variables can be changed and the value of the state variables can be re-initialised at the time that the condition occurs. The pseudo code of a

*when-clause*is the following:where subscript

**represents the value of the variables after execution of the event, subscript**

*a***represents the value of the variables before execution of the event, and is the value of the trigger time of the event. Discrete time variables only change their values at the event times, so they are constant between events.**

*b**Conditional equation clauses*allow changing the active equations in the model. The pseudo code of a conditional equation is:where , and are continuous expressions of the model variables. If the condition is true, then , and if the condition is false, . The tool automatically generates an internal discrete variable

*α*, whose value is one as long as the condition is true and zero as long as the condition is false, and it translates the

*conditional equation*clause into the following equation:

If a *delayed assignment* is triggered, then the values of the *right expression* and the delay are evaluated at the current time and a new *time-event* is scheduled after the delay. Then, the variable is set to the evaluated value of the *expression* when the new *time-event* is triggered.

It is worth noting that after an event is executed, the tool modifies values of the derivatives and the algebraic unknowns ** z** to satisfy the algebraic constraint equations. Whenever an event is triggered, it means stopping the integration to execute the event actions and then restarting it. This process can consume a great deal of CPU time for cases where many events occur.

#### Symbolic manipulation

DAEs are much more difficult to solve than ODEs. OOMS tools manipulate the equations symbolically for two purposes: to free the developer from having to think about how to solve each equation, and to turn the original DAE system into an equivalent but easier to solve DAE, or in some cases into an ODE.

The symbolic manipulation phase carries out a number of different tasks. First, it checks that the number of unknowns is equal to the number of equations and that the equation system is not structurally singular. Second, if there are links between state variables, the number of degrees of freedom is less than the number of states, so the tool automatically reduces the number of states. Third, it attempts to sort the equations into an explicit sequence of calculation, although this is not always possible, sometimes algebraic loops occur. Algebraic loops may be linear or non-linear. Linear algebraic loops are solved by calling a matrix linear solver during simulation. In case of non-linear algebraic loops, the symbolic manipulation attempts to reduce the size of the systems by applying a tearing algorithm and the resulting reduced systems are solved by calling a non-linear solver.

### Simulation environment of OOMS tools

In addition to providing a modelling language, an OOMS tool offers a simulation environment that includes a graphical model editor, an experiment scripting language to define the simulations to be performed, and a graphical tool for plotting and reporting simulation results.

The Modelica-based tools and EcosimPro are examples of object-oriented modelling tools. Modelica is a standard non-proprietary object-oriented modelling language for which several simulation environments are available. Although the idea of a standard language is appealing, it must be noted that portability of models between the different simulation environments is not guaranteed; a Modelica model might work in one environment and fail in another.

### Main disadvantage of OOMS tools

Although OOMS tools make it easier to create simulation applications for different domains, they do have a major disadvantage: they require that the continuous behaviour of the model be described with DAEs. DAEs are a generalization of ODEs and provide a general framework for a simulation tool, but they are not a panacea that can be used for all simulation problems. For instance, transient flow of liquid in pipes is described by two hyperbolic type partial differential equations (PDEs) whose solution consists of the superposition of two waves. This phenomenon is easily captured by the method of characteristics, as it tracks the propagation of the waves. Although it is possible to solve the water hammer PDEs using the method of lines (MOL) by discretizing only in space and then integrating the semi-discrete problem as a system of ODEs or DAEs, the MOC is superior to other methods in the case of water hammer with a constant wave speed (Chaudhry 2014).

## STATE OF THE ART OF WATER HAMMER SOFTWARE

### Commercial water hammer software developed using GPLLs

An extensive review of water hammer theory and practice that includes the description of eleven commercially available water hammer software packages is given in Ghidaoui *et al.* (2005). Of the 11 software packages reviewed, eight use the MOC, two use implicit finite difference methods, and only one uses the wave-plan method. Although not mentioned in the paper, all of them can be assumed to have been developed using GPPLs. AFT Impulse (Applied Flow Technology 2017) is one of these codes and has been used for validation and comparison with the new water hammer code.

### Water hammer software developed using OOMS tools

OOMS tools have also been applied to water hammer simulation. Lindblom (2015) analysed water hammer in a district heating system using the Modelica Fluid library. The purpose was to determine how the tool performed in comparison to the MOC codes when solving water hammer problems. It is concluded that the Modelica library was unsuitable for water hammer simulation because it could only roughly estimate the phenomenon in simple systems.

Koppel *et al.* (2011) documented the implementation and validation of a set of EcosimPro libraries designed to simulate the functioning of a space propulsion system. These libraries are able to simulate water hammer phenomena, particularly the water hammer that occurs during the filling of empty lines.

Although the Modelica and EcosimPro fluid flow libraries used in the two papers cited above are independent developments, both use semi-discrete finite volume methods (FVM) to transform the PDEs that describe the fluid flow into an ODE system. The Modelica library implements a centred differencing scheme whereas the EcosimPro library provides two options: the centred scheme and Roe's upwind scheme.

El-Hefni *et al.* (2011) developed a Modelica library called ThermoSysPro for static and dynamic modelling of power plants. This library provides components to analyse fast dynamic phenomena such as water hammer, and it includes a pipe component that implements the MOC and a simple test case for it. However, because this previous implementation of the MOC does not solve the numerical problems that arise when implementing a discrete solver in an OOMS tool, even the simple test case provided with the library requires a large computational effort.

## METHOD OF CHARACTERISTICS

Classical textbooks on water hammer (Wylie *et al.* 1993; Chaudhry 2014) provide complete derivation of the method of characteristics, but a summary is presented here for the sake of completeness.

### Water hammer equations

The last two terms on the left side of the momentum equation are respectively the quasi-steady friction head loss per unit length and the unsteady friction head loss per unit length. The transient friction term follows the modified instantaneous acceleration model of Vítkovský *et al.* (2006), which was selected among the different unsteady friction models available in the literature because of its simplicity and widespread usage.

*R*and

*S*in Figure 5, Equations (13) and (14) can be solved for the piezometric head and the flow at the new time step in point

*P*, yielding:

To apply the MOC to a pipeline, a rectangular grid such as the one shown in Figure 6 is considered. The pipe length is divided into *N* equal sized reaches (or intervals) of length, *,* and the time is divided into equal time steps, . The edges of the length intervals are designated as ‘nodes’ or ‘sections’. At the initial time, , the head and flow at the nodes are assumed known from a steady state solution. The solution process consists in finding *H* and *Q* for the grid points along then and so on, until the desired time duration is reached. The characteristic lines passing through point *P*, a typical grid point at the next time step, intersect the *x*-axis at *R* and *S*. Although head and flow are only known at the grid points, it is possible to interpolate their values at R and S by using linear or higher order interpolation.

The most obvious selection for the time step is to make the Courant number equal to 1, then points R and S are coincident with the grid points at the previous time step, and the MOC is highly accurate. However, to satisfy for each pipe in a model poses a challenge, because the time step needs to be common but individual pipes feature different properties. This discretization problem can be addressed by artificial adjustments of the wave speed or length, or by interpolation techniques. There are multiple interpolation techniques, linear or non-linear, interpolation in space or in time. Karney & Ghidaoui (1997) derived a discretization algorithm that supports blending of the wave speed adjustment method, the length adjustment method and several types of linear interpolation. Liqhammer has implemented this discretisation algorithm because it constitutes a flexible tool for investigating the importance of discretisation errors in pipeline systems.

### Boundary conditions

At either ends of a pipe, there are two unknowns (, but only one of the compatibility equations is available. For the upstream pipe end in Figure 6, Equation (14) holds along the characteristic, and for the downstream end, Equation (13) holds along the characteristic. Additional relations (equations) between head and discharge must be specified at the pipe ends.

The connecting components located at the pipe ends, which are called ‘boundary elements’, or simply ‘boundaries’, provide the needed additional equations.

Two simple boundary elements – reservoirs and exit valves – are presented next to illustrate how boundary conditions are treated in the MOC.

#### Reservoir

#### Exit valve or orifice

### Column separation and discrete cavitation model

Although the column separation phenomenon is not fully understood, a variety of numerical algorithms have been developed for modelling it in water distribution systems, ranging from single component to two-component two-phase transient solutions. Excellent reviews of the methods can be found in the literature (Simpson & Bergant 1994; Bergant *et al.* 2006).

The discrete vapour cavity model (DVCM) is the most commonly used model for column separation and distributed cavitation. Its advantages are that it is easy to implement and that it reproduces many features of the physical events of column separation in pipelines. Cavities are allowed to grow at any of the computational sections if the pressure becomes lower than the liquid vapour pressure.

## OVERVIEW OF LIQHAMMER CODE

### Programming overview of LiqHammer

The main problem of implementing the MOC using an OOMS tool is to achieve an efficient implementation. The MOC is a discrete method and its interactions with the built-in solvers of the OOMS tools can easily hinder the efficiency of the simulation.

An academic version of the code is available on github, the direction of the repository is given in the Data Availability Section at the end of the paper. Supplementary Material, Appendix A describes some specific parts of the code to help the interested readers on its understanding.

#### Implementation of the calculation of inner sections

The MOC calculation of the inner sections of a single pipe can be implemented in an OOMS tool by means of discrete equations activated by a sample event generator with a period equal to the time step of the MOC.

#### Implementation of the calculation of boundary conditions

In the conventional implementation of the MOC using a GPPL, the boundary conditions are only solved at discrete time steps and the variables are updated discontinuously in time. The same strategy could have been followed when implementing the MOC using an OOMS tool. However, it is preferable to solve the boundary conditions continuously because it prevents numerical discontinuities from being propagated into connected components. A very similar problem occurs in co-simulation, where the discrete update of the coupling variables introduces discontinuities that are propagated to the underlying subsystems and may severely hinder the performance of the local solvers (Andersson 2016).

*t*between the discrete times of the MOC, the characteristic line from the pipe end intersects the

*x*-axis at point

*S*, the piezometric head and flow at

*S*are calculated using linear interpolation between sections 1 and 2, and the compatibility Equation (14) can be applied along the characteristic line

*A–S*.

#### Discretisation and synchronisation of all model pipes

Implementing the MOC with GPPLs requires a common time step used in every pipe to solve the governing equations. Even though implementation of the MOC using an OOMS tool does not require the same time step in all pipes (the tool is theoretically able to manage the events coming from all the model pipes), LiqHammer is designed to consider a common time step to achieve a similar efficiency to that of commercial water hammer software.

If each pipe uses its own integration step, the number of events during the integration increases with the number of pipes in the model. It is known that the variable step solvers typically used for integration in the OOMS tools become very inefficient when multiple events occur (Casella 2015), as each event forces a re-initialisation of the integration process. Using the same time step for all pipes makes the number of discrete events independent of the number of pipes in the model.

From the point of view of programming with an OOMS language, the main difficulty of implementing the discretisation algorithm is the need to loop through all the pipe components. OOMS languages provide information hiding, one component only shares its port variables with other components.

### Users overview of the LiqHammer code

The Pipe implements two optional methods: the MOC and the AWH method (Wylie *et al.* 1993). The AWH method is a variant of the MOC that applies the characteristic equations between the initial and final sections of the pipes avoiding the calculation of the inner nodes. Although it has the disadvantage that its friction calculation is less accurate, it is suitable for the analysis of water distribution systems with thousands of pipes. The MOC option implements transient friction and the DVCM model.

Effective implementation of the MOC and the AWH requires an equal time step across the whole network. Each pipe needs to have at least one reach, and the time step is limited by the Courant number. In case of a system with pipes of very different scales, the simulation becomes extremely expensive because the pipe with the shortest travel time limits the overall time step and the reach length in every other pipe. The component PipeShort avoids the Courant limit on the time step by using an implicit finite difference formulation proposed by Karney (1984) and it can be used in case of pipes with a very short travel time.

A water hammer model is built as a network of pipes interconnected by boundary elements where flow devices and fittings are located. The boundary elements include tanks, collectors, exit valves and inner valves, non-return valves, dead ends, area changes, pumps, air inlet/exit valves, gas accumulators and surge tanks.

A special SolverControl type component has to be included in every LiqHammer model, for several purposes: to calculate the steady state, to discretise the model using the algorithm by Karney & Ghidaoui (1997) and to select the transient solver (MOC or AWH) and its operating parameters.

It is worth noting that thanks to the use of an OOMS language, LiqHammer requires only about 1,500 lines of code. Just for comparison, years ago, one of the authors of this paper participated in the development of a Fortran code for water hammer analysis using the MOC, the code had over 6,000 lines of code without providing any interactive graphical capability.

## EXAMPLE CASES AND COMPARISON TO OTHER WATER HAMMER SOFTWARE

A set of test cases carefully selected from the water hammer literature were defined to verify and validate the various component models in the LiqHammer code. Results for a relevant sub-set of three test cases are presented in this section. These cases were also used to compare the efficiency and accuracy of the new code with other water hammer software.

### Example cases

#### Example case 1: pipe network with instantaneous valve closure

^{3}/s, and the transient is generated by the instantaneous closure of the outlet flow at time . The purpose of the case is to verify the pipe component and the different solver options: MOC or AWH, interpolation scheme and transient or steady friction.

#### Example case 2: pump trip

#### Example case 3: transient cavitation

*et al.*(1993) was selected to verify the implementation of the discrete cavitation model. Figure 16 shows the model, consisting of a pipeline with upward slope in the direction of flow. The inflow is imposed at the upstream boundary and the line is connected to a downstream reservoir whose level is 15 m. The initial inflow in the pipe is 0.4 m

^{3}/s and it is suddenly stopped to initiate the transient at time zero.

### Comparison of the new code to the existing water hammer software

A parametric study of the different methods applied to the example case 1 with different discretisation schemes was performed and the effect of the discretisation on the accuracy and computation time was subsequently assessed.

First, the codes representative of the GPPL and the OOMS approaches were selected. The code selected as representative of the GPPL approach was AFT Impulse, a water hammer software analysis tool that uses the method of characteristics. The reason for this selection is that it is widely accepted by industry and was available to the authors.

Regarding the selection of the OOMS tool, Modelica provides a FLUID library within its set of standard libraries that has been applied to solve some water hammer cases (Lindblom 2015). EcosimPro offers a fluid flow library called FluidFlow1D, within ESPSS (Koppel *et al.* 2011), which is also validated for water hammer analysis. Both libraries apply a staggered finite volume method and can solve two-phase flow. However, it is unfair to use either of these libraries for this comparison because their wider applicability implies a computational overburden in comparison to the new approach – which is limited to classical water hammer. Finally, the PipeliqTran (EAI 2019) library of EcosimPro was the OOMS tool selected for the comparison. This library also applies the staggered finite volume method, but it is specifically designed for solving fluid transients in liquid pipe networks. Supplementary Material, Appendix B presents the staggered finite volume method applied in PipeliqTran.

Before proceeding with the comparison, it is important to note that the tools that implement the MOC, namely LiqHammer and AFT Impulse, automatically discretise the pipes of the model. The user only has to specify the number of reaches in the controlling pipe (i.e., the pipe with the shortest wave travel time) and the rest of the pipes are automatically discretised by the tool. In contrast, with the conventional OOMS tool, PipeliqTran, the user must define the number of control volumes in each pipe. To allow a comparison under similar conditions, the number of control volumes in each pipe of the PipeliqTran model was defined equal to the number of reaches calculated automatically for the same pipe by the tools using the MOC.

The computer used to run the tests is an Intel i5-4570 CPU @ 3.20 GHz, with 8.00 GB RAM. EcosimPro version used is 5.10.2 with the GNU C ++ version 4.7.

. | New Code – LiqHammer OOMS code using MOC . | AFT impulse GPPL code using MOC . | PipeliqTran OOMS using semi-discrete . | ||||||
---|---|---|---|---|---|---|---|---|---|

No. of reaches in controlling pipe . | CPU time (s) . | Max. head at valve (m) . | Error (%) . | CPU time (s) . | Max. head at valve (m) . | Error (%) . | CPU time (s) . | Max. head at valve (m) . | Error (%) . |

4 | 0.14 | 383.4 | 0.58 | 0.09 | 383.5 | 0.66 | 0.34 | 340.7 | −21.74 |

8 | 0.17 | 383.1 | 0.47 | 0.20 | 383.2 | 0.51 | 0.83 | 348.1 | −17.84 |

12 | 0.20 | 383.2 | 0.50 | 0.28 | 383.3 | 0.54 | – | – | – |

16 | 0.25 | 383.1 | 0.47 | 0.39 | 383.2 | 0.50 | – | – | – |

20 | 0.34 | 383.2 | 0.49 | 0.52 | 383.2 | 0.52 | 3.61 | 361.5 | −10.86 |

24 | 0.41 | 383.1 | 0.47 | 0.66 | 383.2 | 0.50 | – | – | – |

28 | 0.50 | 383.2 | 0.48 | 0.82 | 383.2 | 0.51 | – | – | – |

32 | 0.59 | 383.1 | 0.47 | 1.00 | 383.2 | 0.50 | – | – | – |

36 | 0.71 | 381.5 | −0.39 | 1.20 | 381.5 | −0.39 | – | – | – |

40 | 0.83 | 381.6 | −0.33 | 1.46 | 381.6 | −0.32 | 11.60 | 367.8 | −7.57 |

60 | 1.49 | 382.1 | −0.10 | 3.00 | 382.1 | −0.06 | – | – | – |

80 | 2.14 | 382.4 | 0.06 | 4.37 | 382.4 | 0.10 | 36.67 | 372.7 | −5.00 |

100 | 3.03 | 382.5 | 0.13 | 6.29 | 382.6 | 0.18 | 55.18 | 374.2 | −4.19 |

. | New Code – LiqHammer OOMS code using MOC . | AFT impulse GPPL code using MOC . | PipeliqTran OOMS using semi-discrete . | ||||||
---|---|---|---|---|---|---|---|---|---|

No. of reaches in controlling pipe . | CPU time (s) . | Max. head at valve (m) . | Error (%) . | CPU time (s) . | Max. head at valve (m) . | Error (%) . | CPU time (s) . | Max. head at valve (m) . | Error (%) . |

4 | 0.14 | 383.4 | 0.58 | 0.09 | 383.5 | 0.66 | 0.34 | 340.7 | −21.74 |

8 | 0.17 | 383.1 | 0.47 | 0.20 | 383.2 | 0.51 | 0.83 | 348.1 | −17.84 |

12 | 0.20 | 383.2 | 0.50 | 0.28 | 383.3 | 0.54 | – | – | – |

16 | 0.25 | 383.1 | 0.47 | 0.39 | 383.2 | 0.50 | – | – | – |

20 | 0.34 | 383.2 | 0.49 | 0.52 | 383.2 | 0.52 | 3.61 | 361.5 | −10.86 |

24 | 0.41 | 383.1 | 0.47 | 0.66 | 383.2 | 0.50 | – | – | – |

28 | 0.50 | 383.2 | 0.48 | 0.82 | 383.2 | 0.51 | – | – | – |

32 | 0.59 | 383.1 | 0.47 | 1.00 | 383.2 | 0.50 | – | – | – |

36 | 0.71 | 381.5 | −0.39 | 1.20 | 381.5 | −0.39 | – | – | – |

40 | 0.83 | 381.6 | −0.33 | 1.46 | 381.6 | −0.32 | 11.60 | 367.8 | −7.57 |

60 | 1.49 | 382.1 | −0.10 | 3.00 | 382.1 | −0.06 | – | – | – |

80 | 2.14 | 382.4 | 0.06 | 4.37 | 382.4 | 0.10 | 36.67 | 372.7 | −5.00 |

100 | 3.03 | 382.5 | 0.13 | 6.29 | 382.6 | 0.18 | 55.18 | 374.2 | −4.19 |

Figure 18 shows the error of the maximum overhead in the valve versus the number of reaches in the controlling pipe. The new code and the GPPL code offer almost identical accuracy, as expected since they both use the MOC method. In addition, the MOC codes achieve very high accuracy even with a small number of reaches – the error of the maximum overhead in the valve is less than 0.6% with only four reaches in the controlling pipe. The conventional OOMS using ODEs is much less accurate than the MOC and it requires a much finer discretisation to achieve reasonable accuracy – an accuracy better than 5% requires 80 reaches in the controlling pipe.

Figure 19 shows the computation times of the three codes. The new code is the most efficient of them. The computation times of the new code, while lower, have the same order of magnitude as those of the GPPL code. Surprisingly, the new code is slightly faster than the GPPL code. It can also be seen that the computation times of the OOMS approach increase quickly with the number of reaches. This, combined with the need of a fine discretisation to achieve reasonable accuracy, makes the OOMS code inefficient for analysing fast transients. For solving case 1, the OOMS code requires 20 reaches in the controlling pipe to achieve 10% error in the maximum overhead – corresponding to a computation time of 3.61 s – whereas the new code achieves an accuracy of 0.5% with only four reaches and requires a computation time of 0.17 s. Hence, the new code achieves an error twenty times smaller than the OOMS code in less than a twentieth of the time.

## CONCLUSIONS

This paper describes the use of an object-oriented simulation environment to develop an interactive and efficient program for the simulation of hydraulic transients. Interactivity was achieved by using the graphical capabilities of the existing object-oriented simulation environment, and efficiency was obtained by implementing the MOC and the AWH method instead of relying on the built-in ODE and DAE solvers in the OOMS environment.

The adoption of an object-oriented simulation environment provides all the benefits of a generic application framework for development of scientific software (Forde & Stiemer 1989):

*Cheaper development:*LiqHammer was developed much faster using an OOMS tool than using a GPPL. It required only about 1,500 lines of code.*Consistency:*All simulators developed using an OOMS tool share the same simulation environment, namely the language and the graphical user interface. This shortens the learning curve because users can apply the experience gained with a simulator to any other simulation program developed with the same environment.*Transparency:*The LiqHammer source code is easy to understand because the simulation language hides the implementation details.*Flexibility:*Unlike in black box commercial tools, advanced users can modify the code and add new components without any limitation.

Thanks to the use of the MOC, the tool presented here offers much better accuracy and calculation speed than water hammer codes built using an object-oriented simulation environment in a conventional way. The performance of the new tool is similar to or slightly better than the performance of the MOC commercial code selected for comparison, at least for the test cases considered in this study.

The proposed new approach features most of the advantages provided by the OOMS tool: a language that allows the hierarchical definition of models and the definition of the experiments, a graphical interface for drawing the topology of the models and entering the components data and a graphical interface for the analysis of results.

The work presented in this paper was carried out using the object-oriented modelling tool EcosimPro, but it could also have been implemented using the Modelica language. In fact, some initial steps in this direction have already been taken.

## ACKNOWLEDGEMENTS

We thank Empresarios Agrupados International (EAI) who provided the EcosimPro and PipeLiqTran software, and our colleagues in the EcosimPro development team who provided insight and expertise that greatly assisted the research. This work was supported by Ministerio de Ciencia e Innovación of Spain: grant number RTI2018-094665-B-I00 PID2020-112658RB-I00.

## DATA AVAILABILITY STATEMENT

All relevant data are available from an online repository or repositories (https://github.com/RamonVara/REPO_LIQHAMMER_ST).

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

*Methods and Tools for Co-Simulation of Dynamic Systems with the Functional Mock-Up Interface*

*Doctoral Theses in Mathematical Sciences*

*Analysis of Fluid Transients in Large Distribution Networks*

*Ph.D. Thesis*

*Progress in Propulsion Physics*, vol. 2. EDP Sciences, Les Ulis, France, pp. 743–764.

*Modelling of Dynamic Pressure Conditions in District Heating Systems*