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.

  • 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

Graphical Abstract
Graphical Abstract
     
  • pipe cross-sectional area (m2)

  •  
  • valve exit area when fully open (m2)

  •  
  • wave speed of the water hammer waves (m/s)

  •  
  • pipe compressibility constant (s/m2)

  •  
  • quantity for negative characteristics used for 1D MOC solution (s/m2)

  •  
  • quantity for positive characteristics used for 1D MOC solution (s/m2)

  •  
  • 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/s2)

  •  
  • 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/m2)

  •  
  • volume flow (m3/s)

  •  
  • upstream volume flow (m3/s)

  •  
  • downstream volume flow (m3/s)

  •  
  • resistance coefficient (s2/m5)

  •  
  • distance along the pipe (m)

  •  
  • time (s)

  •  
  • vapour cavity volume (m3)

  •  
  • 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/m3)

  •  
  • relative valve opening (unitless)

Subscripts

     
  • ith section

  •  
  • index for pipes connected by inlet end

  •  
  • index for pipes connected by outlet end

  •  
  • last section or number of reaches (depending on the context) of the jth pipe

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.

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.

OOMS languages are declarative, so they can be used to define new port types and new component types. As OOMS languages are probably unknown to most readers, two examples are provided here to illustrate their concepts. The first example is the declaration in EcosimPro language of a hydraulic connection type, the code of which is shown in Figure 1. Port types are defined by a set of variables to be shared between the connected ports and, optionally, by equations relating them. Statements declaring the variables may have prefixes. The SUM prefix for the volume flow variable indicates that the sum of all flows in connected ports must be zero. The EQUAL prefix for the piezometric head indicates that all the heads in connected ports must be equal.
Figure 1

Example of the declaration in EcosimPro of a hydraulic connection type.

Figure 1

Example of the declaration in EcosimPro of a hydraulic connection type.

Close modal
The other example is the declaration in EcosimPro language of a hydraulic component representing a hydraulic exit valve, the code of which is shown in Figure 2. The declaration of a component type includes the declaration of the ports, data and internal variables of the component followed by the continuous equations, which can be algebraic or differential. Optionally, it can also include discrete events and their associated actions, and a definition of an aggregate (a topology), although the example presented does not include them.
Figure 2

Example in EcosimPro of the declaration of a hydraulic exit valve component type.

Figure 2

Example in EcosimPro of the declaration of a hydraulic exit valve component type.

Close modal

Mathematical formulation of OOMS tools

Differential algebraic equations

From a mathematical point of view, the OOMS tools represent the model as a system of DAEs with discontinuities. The general form of a DAE is:
formula
(1)
formula
(2)
where are the derivatives of the state variables are the algebraic variables, are the discrete variables, and 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.

With 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:
formula
(3)
where subscript a represents the value of the variables after execution of the event, subscript b 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.
Conditional equation clauses allow changing the active equations in the model. The pseudo code of a conditional equation is:
formula
(4)
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:
formula
(5)
Delayed assignments can be used as one of the actions within a When-clause. The syntax of a delayed assignment is
formula
(6)

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.

EcosimPro is an OOMS tool with its own modelling language. It also provides a dedicated experiment description language that can set up different simulation scenarios. The scenarios may be simple, such as the calculation of steady states or transient solutions, or complex, such as the solution of optimisation problems. The EcosimPro simulation environment offers a graphical model editor (see Figure 3) that can be used to build models by dragging-and-dropping components from a palette onto a canvas, connecting the ports and entering data on the components; and it also offers a monitoring tool for analysing the results and that can be used for plotting and reporting (see Figure 4).
Figure 3

Model creation window of the EcosimPro simulation environment.

Figure 3

Model creation window of the EcosimPro simulation environment.

Close modal
Figure 4

Monitoring tool of the EcosimPro simulation environment.

Figure 4

Monitoring tool of the EcosimPro simulation environment.

Close modal

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).

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.

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 two equations that describe one-dimensional transient flow of a liquid in a pipe are conservation of mass and conservation of momentum:
formula
(7)
formula
(8)

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.

The MOC transforms the original systems of PDEs into four ODEs:
formula
(9) and (10)
formula
(11) and (12)
Equations (9) and (11) are called compatibility equations and are valid respectively along the characteristic lines defined by Equation (10) and (12). Figure 5 shows the characteristic lines and on the independent variables plane (the plane). The characteristic lines are straight lines of slope .
Figure 5

Characteristic lines in the x–t plane.

Figure 5

Characteristic lines in the x–t plane.

Close modal
By integrating Equation (9) between points R and P, and integrating similarly Equation (11) between S and P yields the following two compatibility equations:
formula
(13)
formula
(14)
where:
formula
formula
formula
formula
is the Karney & McInnis (1990) linearization constant for the friction (. Assuming that the flow and head are known at points 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:
formula
(15)
formula
(16)

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.

For the method to be stable, the time step must fulfil the Courant condition, which requires the Courant number to be less than or equal to 1:
Figure 6

Fixed grid for MOC application.

Figure 6

Fixed grid for MOC application.

Close modal
formula
(17)

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

In the case of a tank or reservoir connected to the upstream end of a pipe (section 1), as shown in Figure 7, there are two unknowns: and . Assuming that the head loss and velocity at the reservoir are negligible, the head at the section 1 of the pipe connected is equal to the piezometric head in the reservoir, Equation (18). The additional required equation is the compatibility Equation (14) along the negative characteristic, .
formula
(18)
formula
(19)
Figure 7

Reservoir.

Exit valve or orifice

In the case of an exit valve or orifice located at the downstream end of a pipes (section ), as depicted in Figure 8, there are also two unknowns: and . The transient flow through the valve is given by Equation (20), and the additional required equation is the compatibility equation along the positive characteristic , Equation (21).
formula
(20)
formula
(21)
Figure 8

Exit valve or orifice.

Figure 8

Exit valve or orifice.

Close modal

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.

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).

It must be noted that the equations representing the boundary condition constraints, such as Equations (18) and (20), are continuous equations which must be fulfilled at any time. The only equations that are discontinuous at the boundaries are the characteristic equations because they are applied at the discrete time steps of the MOC. Turpin (2004) proposed a solution to this problem consisting in the continuous time interpolation of the characteristic equations at the boundaries as shown in Figure 9. At a time 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.
Figure 9

Interpolation of the characteristics at the pipe boundary nodes.

Figure 9

Interpolation of the characteristics at the pipe boundary nodes.

Close modal

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

Figure 10 shows the palette of components of the LiqHammer code. The palette is divided into three groups: pipes, boundary elements and solver control. The Pipes group includes two component types: the Pipe, which is the standard pipe component using the MOC, and the PipeShort, which is designed to represent very short pipes.
Figure 10

Palette of components of the LiqHammer library.

Figure 10

Palette of components of the LiqHammer library.

Close modal

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.

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

First example is a benchmark case from Merilo (1992). Figure 11 shows the model, which includes one reservoir, nine pipes and a variable flow condition representing a valve. The initial condition is a steady state calculated for an outlet discharge of 0.8495 m3/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.
Figure 11

Model of example case 1: pipe network with instantaneous valve closure.

Figure 11

Model of example case 1: pipe network with instantaneous valve closure.

Close modal
Figure 12 compares the results of the new code for the piezometric head at the valve (using the MOC and the steady friction model) with the results calculated by AFT impulse. The results show good agreement and verify the implementation of the MOC for the component types used in the model.
Figure 12

Results of example case 1.

Figure 12

Results of example case 1.

Close modal

Example case 2: pump trip

The second example, shown in Figure 13, is taken from Chaudhry (2014). It represents a pumping system and its purpose is to check the pump component. Initially, the pump is operating at rated condition and the transient is initiated by a pump trip due to a power failure.
Figure 13

Model of example case 2: pump trip.

Figure 13

Model of example case 2: pump trip.

Close modal
Figure 14 shows the results calculated for the fraction of rated speed and the fraction of rated flow during the transient. Figure 15 shows the piezometric head at the pump discharge and at the junction between the two discharge pipes. Both figures also show the predictions reported in Chaudhry (2014) for this case. There is a good agreement between the results of the new code and those from the literature. The code reproduces the successive phases of the transient with great accuracy: the quick pump speed reduction, the flow reversal, the pump running as a turbine and the final pressure rise.
Figure 14

Results of example case 2: dimensionless pump speed and flow.

Figure 14

Results of example case 2: dimensionless pump speed and flow.

Close modal
Figure 15

Results of example case 2: piezometric heads.

Figure 15

Results of example case 2: piezometric heads.

Close modal

Example case 3: transient cavitation

This case taken from Wylie 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 m3/s and it is suddenly stopped to initiate the transient at time zero.
Figure 16

Model of example case 3: transient cavitation.

Figure 16

Model of example case 3: transient cavitation.

Close modal
Figure 17 compares the results of the new code to the results calculated with AFT impulse. The piezometric head at the pipe inlet is practically the same and some minor differences appear only after the third pulse. These minor differences can be attributed to the fact that the DVCM can generate unrealistic pressure pulses (spikes) due to the collapse of multi-cavities.
Figure 17

Results of example case 3.

Figure 17

Results of example case 3.

Close modal

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.

The results obtained for case 1 with the three approaches as a function of the discretisation are summarised in Table 1 and are shown graphically in Figures 18 and 19. The first column of the table shows the number of nodes in the controlling pipe, which is the value that defines the refinement of the discretisation. The subsequent columns show the computation time, the maximum head at the exit valve and the error of the maximum overhead at the valve for each of the three approaches.
Table 1

Comparison of the results of the new code for case number 1

New Code – LiqHammer OOMS code using MOC
AFT impulse GPPL code using MOC
PipeliqTran OOMS using semi-discrete
No. of reaches in controlling pipeCPU 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 (%)
0.14 383.4 0.58 0.09 383.5 0.66 0.34 340.7 −21.74 
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 pipeCPU 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 (%)
0.14 383.4 0.58 0.09 383.5 0.66 0.34 340.7 −21.74 
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

Error of the maximum overhead at the valve for the different codes.

Figure 18

Error of the maximum overhead at the valve for the different codes.

Close modal
Figure 19

Computation time for case 1 using the different codes.

Figure 19

Computation time for case 1 using the different codes.

Close modal
The error of the maximum overhead at the valve is measured using Equation (22), where is the maximum computed head at the valve, is the head of the tank, equal to 191 m, and is a reference maximum head at the valve estimated from a calculation with 400 reaches in the controlling pipe, equal to 382.25 m.
formula
(22)

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.

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.

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.

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

The authors declare there is no conflict.

Andersson
C.
2016
Methods and Tools for Co-Simulation of Dynamic Systems with the Functional Mock-Up Interface
.
Doctoral Theses in Mathematical Sciences
,
Lund University
.
Applied Flow Technology
2017
AFT Impulse 6.0 Waterhammer User Guide
.
Bergant
A.
,
Simpson
A. R.
&
Tijsseling
A. S.
2006
Water hammer with column separation: a historical review
.
Journal of Fluids and Structures
22
(
2
),
135
171
.
Elsevier
.
Casella
F.
2015
Simulation of large-scale models in modelica: State of the art and future perspectives
. In:
Proceedings of the 11th International Modelica Conference
,
September 21–23, 2015
,
Versailles, France
.
Linköping University Electronic Press
, pp.
459
468
.
Chaudhry
M. H.
2014
Applied Hydraulic Transients
.
Springer
,
New York
.
EAI
2015
EcosimPro User Manual
.
EAI
2019
PIPELIQTRAN – Thermo-Hydraulic Transient Simulation
.
Available from: https://www.ecosimpro.com/products/pipeliqtran/ (accessed 22 November 2021)
.
El-Hefni
B.
,
Bouskela
D.
&
Lebreton
G.
2011
Dynamic modelling of a combined cycle power plant with ThermoSysPro
. In:
Proceedings of the 9th Modelica Conference
. pp.
365
375
.
Forde
B. W. R.
&
Stiemer
S. F.
1989
Development of engineering software within a generic application framework
.
Computer-Aided Civil and Infrastructure Engineering
4
(
3
),
205
216
.
Wiley Online Library
.
Ghidaoui
M. S.
,
Zhao
M.
,
McInnis
D. A.
&
Axworthy
D. H.
2005
A review of water hammer theory and practice
.
Applied Mechanics Reviews
58
(
1/6
),
49
.
American Society of Mechanical Engineers
.
Karney
B. W.
1984
Analysis of Fluid Transients in Large Distribution Networks
.
Ph.D. Thesis
,
University of British Columbia
.
Karney
B. W.
&
Ghidaoui
M. S.
1997
Flexible discretization algorithm for fixed-grid MOC in pipelines
.
Journal of Hydraulic Engineering
123
(
11
),
1004
1011
.
American Society of Civil Engineers
.
Karney
B. W.
&
McInnis
D.
1990
Transient analysis of water distribution systems
.
Journal – American Water Works Association
82
(
7
),
62
70
.
Koppel
C.
,
Moral
J.
,
De Rosa
M.
,
Pérez
R.
&
Omaly
P.
2011
Satellite propulsion modeling with Ecosimpro: comparison between simulation and ground tests
. In: Progress in Propulsion Physics, vol. 2. EDP Sciences, Les Ulis, France, pp. 743–764.
Lindblom
A.
2015
Modelling of Dynamic Pressure Conditions in District Heating Systems
.
Master Thesis
,
Chalmers University of Technology
.
Mattsson
S. E.
&
Elmqvist
H.
1997
Modelica-An international effort to design the next generation modeling language
.
IFAC Proceedings Volumes
30
(
4
),
151
155
.
Elsevier
.
Merilo
M.
1992
Water Hammer Prevention, Mitigation, and Accommodation. Volume 4: Review of Analytic Models and Computer Codes. Part 1: Sample Problems and Comparisons. EPRI NP-6766
.
Simpson
A. R.
&
Bergant
A.
1994
Numerical comparison of pipe-column-separation models
.
Journal of Hydraulic Engineering
120
(
3
),
361
377
.
American Society of Civil Engineers
.
Turpin
J. B.
2004
Variable step integration coupled with the method of characteristics solution for water-hammer analysis, a case study
. In
Proceedings of the 52th Jannaf Propulsion Meeting
.
Vítkovský
J. P.
,
Bergant
A.
,
Simpson
A. R.
&
Lambert
M. F.
2006
Systematic evaluation of one-dimensional unsteady friction models in simple pipelines
.
Journal of Hydraulic Engineering
132
(
7
),
696
708
.
American Society of Civil Engineers
.
Wylie
E. B.
,
Streeter
V. L.
&
Suo
L.
1993
Fluid Transients in Systems
.
Prentice Hall
,
Englewood Cliffs, NJ
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC-ND 4.0), which permits copying and redistribution for non-commercial purposes with no derivatives, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Supplementary data