Implementation and Validation of a Free Open Source 1D Water Hammer Code

: This paper presents a free code for calculating 1D hydraulic transients in liquid-ﬁlled piping. The transient of focus is the Water Hammer phenomenon which may arise due to e.g., sudden valve closure, pump start/stop etc. The method of solution of the system of partial differential equations given by the continuity and momentum balance is the Method of Characteristics (MOC). Various friction models ranging from steady-state and quasi steady-state to unsteady friction models including Convolution Based models (CB) as well as an Instantaneous Acceleration Based (IAB) model are implemented. Furthermore, two different models for modelling cavitation/column separation are implemented. Column separation may occur during low pressure pulses if the pressure decreases below the vapour pressure of the ﬂuid. The code implementing the various models are compared to experiments from the literature. All experiments consist of an upstream reservoir, a straight pipe and a downstream valve.


Introduction
Sudden changes to fluid flow in piping may result in the generation of a pressure surge or a chock wave travelling at the speed of sound. The flow changes may be induced by e.g., pump starting/stopping, pipe rupture/leakage, valve opening/closing etc. The phenomenon is referred to as Water Hammer and is most common for liquid flow, due to the low compressibility of water. Probably the most common water hammer is when a downstream valve closes suddenly. When the valve closes, the liquid column upstream is still moving and when the liquid column is slowed down, the kinetic energy is converted to pressure head. The pressure pulse created is noisy and the pressure pulse may have a magnitude which causes pipe rupture or resulting in other equipment damage. Water hammer may be observed in a variety of applications e.g., in power plant steam generation/cooling systems [1,2], district heating system [3], hydro-power [4], irrigation systems [5,6], domestic plumbing [7], pipeline transport of oil and chemicals [8,9] etc. It is important to be able to accurately model the physical phenomenon for correct dimensioning of pipes, adequate choice of valve closing time etc. to prevent damages.
The water hammer phenomenon is fully described by the continuity and momentum equations. Different methods for solving these equations exist, e.g., the Finite Volume Method (FVM) [10],

Water Hammer Equations and Method of Characteristics
The water hammer phenomenon is described by the continuity equation in Equation (1) and the momentum equation in Equation (2) [21,22].
where g is the gravitational acceleration, H is the piezometric head, t is time, A is the cross-sectional area of the pipe, Q is the volumetric flow rate, x is position in the axial direction, J is the friction term, and α is the angle from horizontal. In the equations, it is assumed that the cross sectional area and the wave speed, a, are constant, and that a >> v which means that the convective terms can be neglected. It is also assumed that the flow is one dimensional and that the compressibility and flexibility of the pipe is accounted for in the wave speed as seen in Equation (3).
where ρ is the density of the fluid, K is the bulk modulus of the fluid, E is Young's modulus of the pipe, D is the inner pipe diameter, e is the pipe wall thickness, and c 1 is a coefficient describing the relation to Poisson's ratio, ν.
The method used to solve the set of Partial Differential Equations (PDEs) is the MOC. The MOC transforms the set of PDEs into the four Ordinary Differential Equations (ODEs) seen in Equations (4) and (5).
dQ dt ± gA a dH dt + AJ ∓ g a sin(θ)Q = 0 (4) dx dt = ±a (5) The equations in Equation (4) can be solved along the characteristic lines with the slope determined by Equation (5), as is illustrated in Figure 1. By use of the characteristic lines, point P can be found using point A and B. This is done by the positive characteristic line, C + , corresponding to a positive a and the negative characteristic line, C − , corresponding to a negative a. In this case the characteristic lines are linear, since it is assumed that a is constant. Approximating Equation (4) with finite differences and integrating along the positive characteristic line and along the negative characteristic line yield Equations (6) and (7), where C p and C m are described with Equations (8) and (9), respectively and B = a gA .
C + : H P = C p − BQ P (6) C − : H P = C m + BQ P (7) Inserting Equation (6) into Equation (7) and solving for H P yields: Q P can then be calculated with Equation (11).
Equations (6) to (9) are used in a rectangular grid to simulate the water hammer. The friction term J is evaluated explicitly.
A rectangular grid is chosen over the less computationally intensive diamond grid to enhance the plotting of the pressure distribution. Note that the rectangular grid, which consists of two independent diamond grids, tends to generate small unphysical oscillations, which are visible on the plots.

Steady State
In the previous section the values at point A and B are assumed to be known which means that the steady state of the system has to be determined. In steady state the volumetric flow rate is constant through the pipe. The head through the pipe is calculated with Equation (12).

Friction
The estimation of the water hammer pressure propagation is largely dependent on the friction forces. Traditionally, the friction is modelled as steady or quasi-steady. The steady and quasi-steady friction is able to predict the pressure at the first pressure peak, but as the wave propagates, the dampening of the pressure is not sufficient [23]. For engineering design, the first and thereby largest pressure peak, for single phase flows, is of primary interest. However, there are also specific situations where the secondary pressure peaks are of interest [24]. An example of this is in the cooling system of a nuclear power plant, where a quick restart of the cooling system is necessary. Another example can be quick restart of (off-shore) oil and gas production facilities in order to reduce loss of production etc. If there are sensitive components behind the closed valve, it is necessary to wait for the water hammer to dissipate to a sufficiently low pressure, at which the sensitive components can survive without damage [25].
A method to accurately describe the dampening of the pressure peak is by modelling the unsteady friction. It is complex and more computationally heavy to model unsteady friction compared to modelling steady friction, but unsteady friction modelling is essential for a more precise estimation of the transient behavior [24]. Unsteady friction models fall into three main groups. •

Convolution Based (CB) •
Instantaneous Acceleration Based (IAB) • Extended Irreversible Thermodynamics (EIT) The CB friction models are dependent on the convolution of the history of accelerations with a weighting function. This was first derived by Zielke [26] for laminar flow. Work has been carried out to extend Zielke's CB model to the turbulent regime by Vardy and Brown, assuming a bilinear turbulent viscosity distribution for both smooth pipes [27] and fully rough pipes [28]. Zarzycki et al. [29] derived a more advanced CB friction model for turbulent flow assuming a four-region turbulent viscosity model. The IAB model suggested by Brunone et al. [30] is dependent on the instantaneous acceleration, the convective acceleration, and a damping coefficient. The damping coefficient can be found either empirically as suggested by [30] or as will be used here, the theoretical values derived by Vardy and Brown [27,28]. The determination of the damping coefficient can be done from the initial Reynolds number as a constant or updated on the local Reynolds number, which has been shown to give better results [31]. IAB models with two damping coefficients exist [32,33], but they are not covered herein because there is no analytical way of determining these and therefore they have to be determined empirically.
The EIT derived by Axworthy et al. [34] is a physically well described model, but it has a relaxation time that is usually determined empirically. Storli and Nielsen [35] have derived values based on position dependent coefficients from CB models.
The friction term J consists of a quasi-steady friction term and an unsteady friction term as in Equation (13) [22]. J = J qs + J us (13) The quasi-steady friction term is the skin friction of the pipe and it is modelled using Equation (14).
J qs = f qs Q|Q| 2DA 2 (14) where f qs is the quasi-steady Darcy friction factor which is updated for the local flow. The unsteady friction is modelled with both convolution based friction models, and an instantaneous acceleration based friction model.

Convolution Based Friction Models
The convolution based friction model was developed by Zielke [26] for laminar flow. The unsteady friction is dependent on the convolution of past accelerations with a weighting function as seen in Equation (15).
where µ is the dynamic viscosity of the fluid and W(τ) is the weighting function, which is dependent on the dimensionless time τ, defined as in Equation (16).
All the CB models are defined as in Equation (15) with different weighting functions. The CB models are implemented in a rectangular grid with first order finite difference approximations as in Equation (17) [36].
where i is the spatial index, j is the time index, n ≥ 3 and is the number of time steps from the beginning of the closure of the valve. It should also been noted that in the implementation of the accelerations in the code, the acceleration vector is filled from the bottom of the vector for simplifying the convolution with the weighting function. Equation (17) is evaluated at both the positive and negative characteristics. The weighting function suggested by Zielke is defined as in Equation (18).  [27,28] extended the range of applicability to the turbulent flow regime for Reynolds numbers up to 10 8 with a two-region turbulent viscosity model. In the two-region turbulent viscosity model, it is assumed that the turbulent viscosity varies linearly from the wall to the core region and is constant in the core region. It is also assumed that the viscosity distribution is constant over time. The approximated weighting function derived for a smooth pipe has the form as seen in Equation (19) [27].
Zarzycki et al. [29] utilizes a four-region turbulent viscosity model. The four regions are the viscous sublayer, the buffer layer, the developed turbulent layer and the turbulent core. The weighting function by Zarzycki et al. [29] is approximated as in Equation (20).

Instantaneous Acceleration Based Friction Models
In the instantaneous acceleration based friction model, suggested by Brunone et al. [30], the unsteady friction is dependent on the instantaneous acceleration, the instantaneous convective acceleration and a damping coefficient as in Equation (21), with the sign correction suggested by Vítkovský [37].
where k is Brunone's friction coefficient, which describes the damping of the head, ∂Q ∂t is the local instantaneous acceleration, and ∂Q ∂x is the instantaneous convective acceleration. The implementation of the IAB model has been done explicitly with first order finite differences as in Equation (22) for the positive characteristic and as in Equation (23) for the negative characteristics.
Vardy's shear decay coefficient for laminar flow is defined as in Equation (25) and for turbulent flow, in smooth pipes, as in Equation (26) [37].
Re log(14.3/Re 0.05 ) The analytical solution is based on the CB model suggested by Vardy & Brown, where Vardy's shear coefficient is derived as a limiting value of the unsteady friction coefficient with the special case of constant acceleration.

Boundary Conditions
All the experiments were conducted with a single pipe with a reservoir at the upstream boundary and a valve at the downstream boundary. Therefore, it was necessary to implement the reservoir and valve boundary conditions in the code.

Reservoir
When using a reservoir, it is assumed that the reservoir is large enough, so that elevation changes during operation can be neglected. The head at the reservoir is assumed to be constant, H P = H R , and the flow rate can be isolated from Equation (7).

Valve
When applying a valve, it is important to correctly describe the closure, since it affects both the magnitude and the shape of the pressure peak. It is possible to approximate the valve behaviour using Equation (27) [21].
where τ v is the dimensionless valve closure time, t c is the actual closure time, t is the time, and m is an adjustable constant. If m is set to zero, then the valve is assumed to close instantaneously which yields the maximum pressure peak. The behaviour of τ v is illustrated in Figure 2 for m > 0. When 0 < m < 1, there will be a rapid decrease in flow rate at the beginning and a slow decrease at the end of the closure. For m = 1 there will be a linear fall in flow rate during the closure. When m > 1, there will a slow decrease in flow rate at the beginning and a rapid decrease at the end of the closure. The flow rate can be calculated with Equation (28).

Column Separation and Cavitation
If the pressure during a low pressure period for a liquid-filled pipe reaches the vapour pressure, the water starts to evaporate. This evaporation will form vapour bubbles and cavities. A subsequent increase in pressure will cause the bubbles to collapse. This phenomenon often referred to as cavitation or column separation has been implemented in the code by two different models-the Discrete Vapour Cavity Model (DVCM) [15,21] and the Discrete Gas Cavity Model (DGCM) [38].
Both the DVCM and the DGCM models assume that the wave speed is not affected by the amount of free gas. It is also assumed that all the free gas in each reach is present in a single pocket at the node causing the wave speed to be constant and identical to the wave speed for a single phase system between each node. The gas pocket at the node is modelled with Equation (29) [21].
where Q out is the volumetric flow rate exiting the node and Q in is the volumetric flow rate entering the node. To determine the volume of the gas pocket, Equation (29) can be integrated from time t − 2∆t to time t, and the result is given in Equation (30) [21].
where V g is the volume of the gas pocket, Q is the volumetric flow rate, the subscript P indicate points at time t, P0 indicate points at time t − 2∆t, u indicate volumetric flow rate entering the node, and ψ is a weighting factor. The volumetric flow rates are illustrated in Figure 3. It can be seen that there is a modification to the negative characteristics equation, C − , since it goes from Q u,B and not Q B . Therefore the Q B in Equation (9) is replaced with Q u,B and Q u,P is calculated by Equation (31). The weighting factor ψ controls the weight of the values at t and at t − 2∆t. With ψ = 0 only values at t − 2∆t are used and with ψ = 1 only values at t are used. It is generally recommended to keep ψ above 0.5 as an over reliance of old values have shown to give excessive numerical oscillation. A value close to 0.5 is expected to give the most accurate results, but at this low value some numerical oscillations may still be present. The numerical oscillation is primarily present when gas cavity sizes are small, i.e., in the high pressure zones of the water hammer event. To decrease the numerical oscillations, the value of ψ can be increased towards unity. However, this will cause more spreading of rarefaction waves which can result in increased attenuation. Setting ψ = 1 will remove all numerical oscillations, but will introduce the largest amount of attenuation. Therefore it is recommended to choose a value of ψ as close to 0.5 as possible while experiencing minimal oscillation [21].

Discrete Vapour Cavity Model
DVCM is the simplest of the two phase models included in the program, and it can be used for many flow conditions [21]. However, as will be seen later, DVCM has its limitations when the void fraction becomes too high. Therefore it is recommended to only use DVCM when the void fraction is below 10% [39].
It is assumed that there is no free gas in the system, and that at steady state, and when the pressure is above the vaporization pressure, there is no vapour present in the system, V cav = 0. The flow can therefore be treated as a single phase system, where the volumetric flow rates can be set equal to each other, Q u,P = Q P , because the difference between them describes the increase in vapour volume for each time step. When the pressure becomes lower than or equal to the vaporization pressure, the nodes are treated as boundary nodes, with a fixed pressure as described in Equation (32) [21].
When the head in the node is known, it is possible to calculate the volumetric flow rates and the vapour cavity size with Equation (11), (30) and (31) respectively, where V g,P and V g,P0 are replaced with V cav,P and V cav,P0 .

DVCM-Steady State
As described earlier, DVCM assumes that flow at steady state can be treated as a single phase system, with V cav = 0, and it is therefore calculated as described in Section 2.2.

DVCM-Boundary Conditions
The upstream reservoir is calculated as in Section 2.6.1, with V cav = 0, because it is assumed that the pressure in the reservoir is always above the vaporization pressure.
For the downstream valve, the outlet volumetric flow rate, Q P , is described as in Section 2.6.2, while the nodes are treated in a similar fashion as in Section 2.7.1. The difference between the calculation method for the interior nodes, see Figure A2 in Appendix A.1, and the valve boundary, is that Equation (10) is replaced with Q u,P = Q P and Equation (33) in the flowchart, while Equation (11) is replaced with Equation (28) in the flowchart in Appendix A.1.

Discrete Gas Cavity Model
DGCM is the more complex model, of the two phase models included in the program, and it can be used for many flow conditions.
In DGCM, it is assumed that there is always a small amount of free gas present in the system, and because there is always a small amount of free gas present, it is not possible to calculate the flow as in DVCM. A new expression for the head, which takes into account the effect of the gas cavity, has to be implemented. This is done via Equation (30), where Equation (31) is used to describe Q u,P and Equation (28) is used to describe Q P . This results in Equation (30) with two unknown parameters, H P and V g,P . V g,P can be described with the ideal gas law, as seen in Equation (34), where it is assumed that the mass of free gas, M g , is constant.
R g is the gas constant, T is the temperature, P g is the absolute partial pressure of the free gas, and P g,0 is the absolute partial pressure for the initial void fraction α 0 (i.e., at steady state). V g,P is isolated from Equation (34) and a simplified form can be seen in Equation (35), where P g is described with Equation (36).
C 3 is a constant which is calculated with Equation (37).
With the expression for V g,P in Equation (35), Equation (30) will have the form seen in Equation (38), remembering that Q u,P and Q P are described with Equations (31) and (28) respectively.
where B 2 , C 4 , B v , and B 1 are defined as in Equation (40).
Equation (39) can be solved as a quadratic equation, where x is the variable. The result can be seen in Equation (41), where the two first expressions are isolated from the roots of the quadratic equation. The third and fourth expressions are linearised versions of the two first expressions, which can produce inaccurate results in extreme conditions of high pressure and very low void fraction, or at very low pressure and high void fraction, where |B B | = |C 4 /B 2 1 | << 1. The fifth expression is for the case when B 1 = 0.
With the head described as in Equation (41), the volumetric flow rates and the gas cavity size are calculated with Equations (31), (11), and (30) respectively.

DGCM-Steady State
For the steady state conditions, DGCM uses the same calculation method as for single phase flow, Section 2.2, with the addition of the calculation of the gas cavity size. At steady state Equation (34) can be simplified to Equation (42) because P g,0 = P g .

DGCM-Boundary Conditions
For DGCM, the upstream reservoir uses the same calculation method as for single phase flow, Section 2.6.1, with the addition of the calculation of the gas cavity size. The gas cavity size is calculated with Equation (35), where H P is replaced with H r .
For the downstream valve, the same method for calculating the outlet volumetric flow rate, Q P , as in Section 2.6.2 is used. The expression for the head is derived in a similar fashion as to Section 2.7.2, but this time an expression for Q P is not needed in Equation (30), as it is now described with Equation (28). This gives the expression in Equation (43), which can be rearranged in Equation (44).
where B 2 , C 4 , B v , and B 1 are defined as in Equation (45).
It can be seen that Equations (44) and (39) has the same form and therefore they are solved in the same form, and give the same results, Equation (41). It is important to note that the expressions for the coefficients in Equation (45) are different from the coefficients in Equation (40).
The volumetric flow rates and the gas cavity size are calculated with Equations (31), (28) and (30), respectively.

Experimental
Previous articles have compared CB and IAB friction models to a single experiment, where a friction model is recommended for a single experiment [40,41]. In this article, three CB friction models and one IAB friction model will be compared to three different experiments to establish how the different friction models behave under different experimental conditions. As for the implementation of different friction models, the implemented cavitation models are compared to experiments from the literature.
The experiments chosen consist of a single pipe between a reservoir/pressure tank at the upstream end and a fast closing valve at the downstream end as seen in Figure 4. The different friction models are compared with experimental results found in literature. The pressure is evaluated at the valve position. The selected experiments vary in pipe material, diameter, length, thickness as well as flow velocity in an attempt to see how these differences affect, prediction of the water hammer phenomenon. An overview of the different parameters is found in Table 1. It is attempted to cover a wide variety of flow conditions to test the implemented models. The three chosen experiments for validating the different implemeted friction models have a Reynolds number of 52052, 15843, and 8453, respectively. The experiments used for comparing the different cavitation models have Reynolds number of 6605, 9894 and 30823, respectively.

Experiment 1
The experimental study conducted by Traudt et al. [42] consist of a 7.7 m stainless steel pipe (grade 1.4541) with a diameter of 19 mm and a wall thickness of 1.5 mm. Young's modulus is set to 200 GPa, the pipe roughness to ε = 2 × 10 −6 m, and Poisson's ratio to ν = 0.28. The test section has a high-pressure tank at the upstream end and a fast acting valve, with a closure time of 18 ms, at the downstream end. The piping also has a 1 • upwards slope. The initial flow velocity is 2.75 m/s with a head at the reservoir of 433 m. For the simulation, m is set to 3.8, which was the best fit for the valve closure, as τ v is not reported by Traudt et al. The fluid in the study is water at room temperature.

Experiment 2
The experimental study conducted by Adamkowski and Lewandowski [40] consists of a 98 m copper pipe in a spiral coil with an inner pipe diameter of 16 mm and a wall thickness of 1.0 mm. Young's modulus and Poisson's ratio are reported in the article at 120 GPa and 0.35 GPa, respectively. The roughness was not stated by [40]. For the present investigation, a roughness of 1.5 × 10 −6 m is used corresponding to a typical value of copper. The test section has a high-pressure tank at the upstream end and a quick-closing spring driven ball valve, with a closure time of 3 ms, at the downstream end. The test section has a reported maximum inclination angle of less than 0.5 • . Because of the uncertainty of the inclination and its small size it is neglected in the simulations. The initial flow velocity is 0.94 m/s with a head at the reservoir of 128 m. For the simulation, m is chosen as 0.1, which was the best fit for the valve closure as τ v is not reported by Adamkowski et al. The fluid in the study is water where the density, viscosity, and bulk modulus was given at 1000 kg/m 3 , 0.9493 × 10 −3 kg/m· s, and 2.1 GPa, respectively.

Experiment 3 and 4
The experimental study conducted by Soares et al. [41] consists of a 15.22 m straight copper pipe with an inner diameter of 20 mm and a wall thickness of 1.0 mm. The test section has a hydropneumatic tank at the upstream end and a pneumatically actuated quarter turn ball valve at the downstream end. The closure time for the valve is not described and therefore it is approximated by the time it takes to reach maximum pressure from steady state pressure. The closure time is estimated at 18 ms. Two experiments were conducted by Soares et al. with an initial velocity of 0.423 m/s (Experiment 3) and 0.497 m/s (Experiment 4). The lowest velocity resulted in a single phase water hammer, whereas the largest velocity resulted in a two phase water hammer. The head at the reservoir is 46 m for the two experiments. For the simulation, m is chosen as 3, which was the best fit for the valve closure as τ v were not described by Soares et al. The material properties for the copper pipes were not given in the study, and the properties from experiment 2 were used. The fluid used in the study is water and the experiment is conducted at 20 • C.

Experiment 5 and 6
The experimental study conducted by Bergant et al. [43] consists of a 37.22 m straight copper pipe with an inner diameter of 22.1 mm and a wall thickness of 1.63 mm. The test section has a pressurized tank at the upstream and the downstream end of the pipe. To generate the transient, a fast closing valve is placed at the downstream end of the pipe. The pipe has an inclination of 3.12 • . Two experiments have been conducted with an initial velocity of 0.30 m/s (experiment 5) and 1.40 m/s (experiment 6) which resulted in column separation. The closure time for both experiments in 0.009 s and the head at the upstream pressurized tank is 22 m. The material parameters were not given by [43], therefore the same parameters used for the [41] experiments is used for experiment 5 and 6 as the pipes consists of the same material. For the simulation, m is set to 5, which was the best fit for the valve closure, as τ v is not reported by [43]. Water is the working fluid at room temperature.

Results and Discussion
In this section the results for the single phase simulations of experiments 1, 2 and 3 and the column separation simulations of experiment 4, 5 and 6 will be presented and discussed. In addition to this, the code has been benchmarked against computational fluid dynamics (CFD) simulations. The results are presented elsewhere [44]. Good agreement was found both for with and without column separation when unsteady friction modelling was applied in the MOC code.

Single Phase
A grid independence study showed no difference for the steady and quasi-steady friction models. For the unsteady friction models the head varied from 12 to 24 reaches, but hereafter the change was insignificant. Therefore all single phase simulations are performed with 24 reaches. Figure 5 and Table 2 show the results of the different friction models for experiment 1. The mean oscillation frequency of the experiment is found at 43.14 Hz which is lower than the frequencies calculated with the Steady and Quasi-Steady friction models with the fastest frequency of 45.14 Hz, followed by Brunone with 44.78 Hz and the CB friction models with a frequency of 44.67 Hz. This causes a phase offset between the experimental and simulated data, which is evident in Figure 5 after the third high pressure peak. At the end of the data the phase offset is almost half an oscillation.
It can be seen that all the friction models have the same tendency at the first peak and are able to accurately estimate the head of the peak with only a slight overestimation. At the third pressure peak, the CB friction models accurately model the friction with an overestimation of the head in the range of 0.43-0.50%, whereas the other friction models start to underestimate the friction. At the tenth pressure peak, none of the friction models accurately describe the dampening of the head with the Steady and Quasi-Steady friction models overestimating the head with 24.64% and 24.23%, the Brunone friction model with 15.33% and the CB models with 10.33-10.52%. From Table 2, it can be seen that the CB based models gives the best prediction of the experimental results.  The experimental results [42] indicate some double peaks, which are not predicted by any of the friction models. This is believed to be caused by harmonics of the experimental setup. As such Traudt et al. [42] report that the second harmonic frequency of the experimental setup of 132 Hz is interfering with the frequency of the water hammer which is 43 Hz. Since none of the models accurately describe the experimental data, a sensitivity analysis on the parameters, Young's modulus, Poisson's ratio and the temperature of the water were conducted to see, whether this could explain the difference. Both Young's modulus and Poisson's ratio are varied ±20% and the temperature of the water ±10 • C. Decreasing the Young's modulus and the temperature and increasing the Poisson's ratio individually improved the results slightly compared to the experiment, but not enough to explain the difference. Investigating combined effects also showed an improvement compared to the experimental data, but again it could not explain the difference. Therefore, there must be an effect in the experiment that causes extra damping on the head and a slower oscillation frequency that is not taken into account in the friction models.
In Table 2, an average calculation time, for five identical simulations, for each friction model can be seen, together with a normalized calculation time, where the reference time is the average calculation time for the steady friction model. It can be seen that as the complexity of the friction model increases, so does the calculation time, where a significant increase is seen when going from Brunone's unsteady friction model to the CB models. The reason for this increase is because of the convolution in the CB models, which uses all of the calculated volumetric flow rates in the used nodes, while Brunone only uses the volumetric flow rate from the two previous time steps. Figure 6 and Table 3 show the results of the different friction models for experiment 2.
The oscillation frequency for all the models is close to the experimental, see Table 3. Looking overall on the first peak, the steady, quasi-steady and Brunone coincide best with the data, and all three CB friction models seem to overestimate the head slightly, but it is close to the highest peak in the experimental data with a difference of −0.01% for Vardy and Brown, 0.63% for Zielke and 0.91% for Zarzycki. The head at the valve in steady state is lower for the experiment compared to the simulations, which could be because of an underestimation of the steady state friction. Initially, the friction models suggested by Zielke and Zarzycki render the most accurate representation of the dampening of the head and at the third pressure peak, only underestimating the head by 0.05% and 0.58%, respectively. This is followed by the friction model by Brunone with an overestimation of 1.01%, the Quasi-Steady model with 1.83%, the Steady model with 2.15%, and the model by Vardy and Brown with 2.35%. At the tenth pressure peak, the picture is different with the model by Vardy and Brown giving the most accurate representation with an overestimation of the head by 1.17% followed by the model by Brunone with an overestimation of 2.31%, the model by Zielke with an underestimation of 3.86%, the model by Zarzycki with an underestimation of 5.31%, the Quasi-Steady model with an overestimation of 7.64%, and the Steady model with an overestimation of 9.81%. The behaviour of the steady, quasi-steady and Brunone models is not reminiscent of the experimental results. The CB based models behaviour is reminiscent of the experimental data with the Vardy & Brown model being the closest to the experimental. The difference in the CB friction models is attributed to the difference in the weighting functions and their dependency on the Reynolds number. Overall, it is the model by Vardy & Brown that provides the best representation.  In Table 3, the average calculation time is summarised. It can be seen, that as the complexity of the friction model increases, the calculation time increases, which was also seen for experiment 1. Figure 7 and Table 4 show the results of the different friction models for experiment 3. The oscillation frequency of the models, ranging from 20.74 to 20.94 Hz, is close to the experimental one, at 20.69 Hz, therefore no offset between the experimental data and the simulations is expected, and it is also evident in Figure 7. The head at steady state is larger for the simulations than the experiment, which could be due to an underestimation of the steady state friction.  All the friction models slightly overestimate the first pressure peak with steady and quasi-steady giving the best estimation with an overestimation of 2.29% and 2.28%, respectively, Brunone being comparable with an overestimation of 2.31% and the CB friction models with an overestimation ranging from 2.87% to 2.89%. Part of the overestimation of the models is likely caused by the overestimation of the steady state head. The model by Brunone represents the dampening of the head most accurately by overestimating the head with 1.67% at the third pressure peak and with 0.59% at the tenth pressure peak, followed by the CB friction models overestimating the third pressure peak in the range 2.05-2.19% and for the tenth pressure peak 1.70-2.28%, the Steady and Quasi-Steady models overestimating the third peak with 5.58% and 5.55%, and for the tenth peak with 17.29% and 17.13%.
In Table 4, the average calculation time is summarised. The same trend as for experiment 1 and 2 is observed.
For all the experiments, it can be seen that the steady and quasi-steady friction models accurately estimate the first peak of the pressure wave, but do not accurately describe the damping of the pressure wave. There is almost no difference in the head of steady and quasi-steady, which means that there is little to no gain from calculating the friction factor for each node and time step which is done in the quasi-steady model. If a precise estimation of the wave propagation is wanted, an unsteady friction model shall be used. A general behaviour for the experiments with copper pipes (experiments 2 and 3) is that the friction in the pipe in steady state is underestimated. This can be due to some of the material properties chosen for the simulations not exactly matching the real properties.
Based on the four experiments, the Vardy & Brown friction model is recommended as a generic choice. The friction model has a good estimation of the wave behaviour and propagation. The head of the friction model is in neither of the experiments underestimated, as e.g., Zielke and Zarzycki in experiment 2.

Column Separation
For both DVCM and DGCM for all three experiments, a grid test was made and it showed that a grid of 24 to 48 reaches was sufficient. Compared to single phase model, an additional parameter (the weighing factor, φ) had to be tuned for the two phase models. It is recommended that φ is set to 0.5 and it was found that values close to this was adequate. It was also found that the best results were obtained for both column separation models in combination with the friction model suggested by Vardy & Brown. Further information of the grid test, choice of φ, and results obtained by the other friction models are found in [44]. Figure 8 and Table 5 show the results obtained with DVCM and DGCM for experiment 4.  The oscillation frequency of both DVCM and DGCM are close to the one obtained from the experimental data of 20.24 Hz. In the first high pressure zone it can be seen that DVCM and DGCM give almost identical results both overestimating the pressure by 2.93%. In the second high pressure zone the large pressure peak is caused by the implosion of bubbles. Both models give an overestimation of the second high pressure peak with DVCM giving an error of 6.28%. In the third pressure zone neither of the models are able to model the largest pressure peak with DVCM underestimating the pressure peak by 30.07% and DGCM by 21.91%. However, both models have a maximum pressure exceeding this pressure in the previous pressure peak, and is therefore still conservative although slightly inaccurate. On the tenth peak DVCM gives the best results overestimating the head by 10.17% compared to DGCM by 17.53%. Overall the DVCM is the best model for experiment 4.
In Table 5, it can be seen that the calculation time for DVCM and DGCM is very similar, and hence the added complexity of DGCM does not extend the calculation time significantly. Figure 9 and Table 6 show the results obtained with DVCM and DGCM for experiment 5.
The oscillation frequency is larger for both DVCM, 8.70 Hz, and DGCM, 8.46 Hz, than the one obtained by the experimental data, 8.15 Hz. Thus, the models are expected to oscillate faster than the experiment, which is also clear for DVCM. For DGCM it seems that there is almost no difference between the oscillation of the model and the experiment. This means that the oscillation of the pressure wave is accurate but that the timing of the pressure peaks are off.  In the first high pressure zone both models give similar results overestimating the pressure by 4.83% or 4.84%. For the second pressure peak both models give a good approximation of the pressure with DGCM giving the best approximation with an overestimation of 4.53%. In the third high pressure zone there is two high pressure peaks; one in the beginning and one in the end. Neither of the models are able to model the pressure peak in the beginning of the third pressure zone, but both are able to model the pressure peak in the end. This pressure peak is most accurately modelled by DGCM with an overestimation of 7.61% compared to an underestimation of 13.80% with DVCM. In the seventh pressure zone DVCM gives the best results with an overestimation of the pressure by 11.96% compared to 19.57% obtained with DGCM. Overall DGCM gives the best result as it accurately model the oscillation of the pressure wave and provides the best results when cavitation is occurring.
In Table 6, it can be seen that the calculation time for DVCM and DGCM is very similar, as was the case with experiment 5, and hence the complexity of DGCM does not extend the calculation time significantly. Figure 10 and Table 7 show the results obtained with DGCM for experiment 6.  In this experiment only the DGCM gives satisfactory results. The DVCM model has a way too fast oscillation and dampening of the pressure, giving unrealistic results. Therefore, the results obtained from DVCM have been omitted. This failure to obtain realistic results is attributed to the large void fractions obtained, which is known to cause problems for the DVCM.
Comparing the oscillation frequency of the DGCM with the experiment it is clear that they are almost identical at 2.92 Hz and 2.90 Hz respectively. Therefore almost no offset between the experimental data and the simulated data are expected, as is also evident in Figure 10. The DGCM give accurate results for the first three pressure zones never overestimating the pressure by more than 5.83%. However, at the fourth pressure zone the simulated data starts to exhibit some oscillations causing an overestimation of 20.09%.
In Table 7, the calculation time for DGCM is shown. The reason for the calculation time being lower than experiment 5, while still with a higher flow time, is due to experiment 5 using 48 reaches, while experiment six uses only 24 reaches. The overall best model for cavitation/column separation is considered to be the DGCM. It is more robust and provides accurate results of the pressure while never underestimating the maximum pressure. DGCM is recommended as a generic choice.

Conclusions
In this paper an implementation of a code for simulation of pressure transients in liquid-filled piping both with and without cavitation/column separation has been presented. The solution scheme implemented for simulating water hammer is the MOC.
For single phase/ liquid-filled pipe water hammer, six friction models are implemented and tested against three different experiments found in literature. The friction models consist of a steady, quasi-steady, and four unsteady models. The unsteady friction models consist of three CB models (Zielke, Vardy and Brown, and Zarzycki) and one IAB model (Brunone). For the experiments included in the present study the Vardy and Brown model seems to provide the best results.
Two models for taking the effects of cavitation/column separation into account have been implemented, the DVCM and DGCM models. The two models are compared to three different experiments. For all experimental comparisons the Vardy and Brown CB friction model has been applied. For the included experiments the DGCM gives the best results overall, although the DVCM provides a slightly better fit for a single experiment.
Generally, it is found that employing unsteady friction models the presented code is adequately able to simulate the chosen experiments.
The authors provide the full verbatim source code to the presented MOC implementation along with the present paper. In doing so, the authors hope that the code can find good use for others. Acknowledgments: Language secretary Susanne Tolstrup, Field Development, Ramboll Energy, has kindly provided her assistance for proofreading this manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following common symbols are used in this manuscript:

. MOC Calculation Flow
In the Octave/MATLAB program all equations are implemented explicitly and a flowchart of the program is in Figure A1. The program starts with importing user defined parameters and settings. The steady state conditions is calculated before the closure of the valve. Then the spatial index, i, and the time index, j, are set to 2. All the interior nodes are calculated, i.e., from i = 2 to i = n − 1 for time j = 2. When i = n − 1 the boundary conditions, i.e., the flow conditions at the reservoir and the valve. Then the accelerations are calculated if an unsteady friction model is used. If the time step j < n t , the time step is updated to j = j + 1 and i = 2 or if j = n t the simulation is finished and the data can be outputted.
A flow chart of the calculation method for the interior nodes for the DVCM can be seen in Figure A2, where the block with V cav (i, j − 1) > 0 investigates the presence of a vapour cavity in the previous time step. If a vapour cavity was present, it is assumed that the current time step should be treated as a pressure boundary. If the node is calculated as a pressure boundary, the calculations are followed by a check for whether the vapour cavity is less than or equal to zero. If this is true, it is assumed that the vapour cavity has condensed, but the pressure has not risen above the vapourization pressure.

Appendix A.2. Code Organisation
To make the program easily adjustable, the code has been divided into a input-file, master-file, solver-file, function-files, boundary condition-files and output-file. All the files are executed in the master-file which is the file that has to be "run". All the user specified values are entered into the input-file, which is also where the solver type, boundary conditions and friction models are chosen. The following boundary conditions have been implemented into the program: Reservoir, instantaneous valve closure and transient valve closure. If additional boundary conditions are wanted, it is easy to define in the input-file, when the boundary condition function-file has been constructed. Additional interior nodes have to be implemented into the solver. The program solves for the setup seen in Figure 4 where the system parameters are specified in the input-file.
The file structure and organisation of code is provided below.