Numerical simulation of water hammer in low pressurized pipe : comparison of SimHydraulics and Lax-Wendroff method with experiment

Article describes simulation of unsteady flow during water hammer with two programs, which use different numerical approaches to solve ordinary one dimensional differential equations describing the dynamics of hydraulic elements and pipes. First one is Matlab-Simulink-SimHydraulics, which is a commercial software developed to solve the dynamics of general hydraulic systems. It defines them with block elements. The other software is called HYDRA and it is based on the Lax-Wendroff numerical method, which serves as a tool to solve the momentum and continuity equations. This program was developed in Matlab by Brno University of Technology. Experimental measurements were performed on a simple test rig, which consists of an elastic pipe with strong damping connecting two reservoirs. Water hammer is induced with fast closing the valve. Physical properties of liquid and pipe elasticity parameters were considered in both simulations, which are in very good agreement and differences in comparison with experimental data are minimal.


Introduction
Unsteady flow in pipeline systems is well known phenomenon.It is possible to describe it by momentum and continuity equations, which can be simplified to solve one dimensional problem.One can find detailed derivation in [1].
There are plenty of numerical methods, which allow solving this hyperbolic problem.Method of characteristic is used in [2] and it is probably the most popular method in fluid mechanics, Lax-Wendroff is chosen in [3], SIM-SEN software [4] uses approach by Runge-Kutta and many other methods are derived in [5].
Behaviour of pipeline system is defined by many factors and sound speed is one of them, probably, one of the most important, see equation 1.
As hydraulic systems does not consist only of fluid, it is necessary to consider pipe properties (Young's modulus, diameter and thickness of the wall), which reduce bulk modulus and so sound speed.When static pressure is low enough it is good idea to consider also influence of undissolved gas in fluid.Low enough could mean under 0.5 MPa, but it depends on content of gas.Influence of air on sound speed in water was measured by [6] and [7] and analytical derivation is published in [8].This paper deals with water hammer in simple pipe with strong damping.It is not possible to neglect influence of undissolved air in water because of low static pressure.Usually, software for computation of water hammer offers some cavitation model, but influence of bubbles of air is not included so numerical simulation do not agree with real a e-mail: daniel.himr@vsb.czresult.Generally, this is problem of programs, which are based on method of characteristic.In this case, two pieces of commercial software were used in this case: SimHydraulics, product of Mathworks Company, uses ode15s numerical method and HYDRA, product of Brno University of Technology, is based on Lax-Wendroff method.Both of them allow considering pipe properties and content of free air in water.

Fundamental Equations
Momentum and continuity equations are sufficient for description of unsteady behaviour of liquid.When simplifying assumptions are introduced, equations have 1D shape respectively (2) and (3).Main assumptions are: -Diameter of pipe does not change in axial direction -Velocity of flow is perpendicular to cross-section of pipe -Pressure is the same in the whole cross-section -Density of fluid does not change Equation (2) includes, beside friction loss, also influence of gravity (right side).Continuity equation (3) contains general function f (p, t), which serves for description of pipe wall pulsations (diameter change as function of pressure).This function can be variable with time.When pipe material is ideally rigid f (p, t) = 0. Typically, when water flow is solved, density is considered to be 1000 kg/m 3 and bulk modulus around 2.1 GPa, but, in fact, water or any other liquid contains some amount of undissolved gas, which can play important role.When pressure drops down, volume of gas bubbles increases and slightly lowers total density of the mixture (fluid+gas) and, at the same time, dramatically influences bulk modulus.Figure 1 shows dependence of density (dashed line) and bulk modulus (solid line) on the pressure.Plotted curves are valid for mixture, which consists of water and air.As influence of pipe wall elasticity is included, maximal bulk modulus in not 2.1 GPa, but lower.Detailed derivation is presented in [8].It leads to definition of sound speed in fluid as function of pressure, see figure 2 and equation 1.Maximal velocity of shock wave is approximately 300 m/s, because of high elasticity of pipe.It is also possible to use 2D or 3D mathematical model to solve flow in pipeline.First approach is described e. g. in [9] and Fluent could be named as representative of 3D way, but this is not suitable when large pipeline system has to be solved.
Following part of the paper is more focused on solutions by HYDRA and SimHydraulics.

SimHydraulics solution
Matlab-Simulink-SimHydraulics is developed to solve hydraulic systems, which are modelled as block schemes.It includes models of different elements like valve, pump, surge tank and so on.Of course, viscoelastic pipe is listed.
It is called segmented pipe and each segment consists of resistance (represents friction loss), fluid inertia and air chamber (represents pipe elasticity and damping).The more segments the more accurate result.
Pipe elasticity is defined via "Static pressure-diameter coefficient" that means diameter change with the pressure.It has unit m/Pa and it is possible to derive relationship between this coefficient and Young's modulus of the pipe material, see equation (4).Formula can be used when pipe wall is thin and material is subject to Hook's law or if diameter change is low enough.
Damping is defined by "Viscoelastic process time constant" that describes velocity of diameter change.When pressure in the pipe increases, diameter changes to a new value D N , but it needs certain time to reach it.Viscoelastic process time constant τ is explained in the figure 3.If Voight-Kelvin model of solid is considered, constant can be expressed as function of damping and elasticity, see equation 5.
Fluid is defined by density, kinematic viscosity, bulk modulus at atmospheric pressure without gas and volume fraction of trapped air at atmospheric pressure.
To solve equations ( 2) and (3) SimHydraulics offers ode15s method, which has adaptive time step and is good for solving stiff problems.

HYDRA solution
HYDRA software uses Lax-Wendroff method, which is based on the Taylor series up to 2nd derivative.Formula for expansion of general function s(x, t) is written in equation 6.  3) are solved in discrete points with constant space step along the pipe.
To define fluid properties, one needs density and kinematic viscosity of liquid and theoretical sound speed (it means 1500 m/s in case of water).Content of gas is described by temperature, mass fraction, gas constant and adiabatic exponent.Relationship between mass fraction and content of air according to SimHydraulics can described by equation (7), which is valid for little volume of free air.

Experiment
SimHydraulics and Hydra were used to simulate water hammer in a simple gravitational pipeline in laboratory.The circuit consists of two tanks, which are connected by polyethylen pipe.Downstream end of the pipe can be closed with a valve.Overflow ensures that head is constant (see fig. 4).Important properties ate listed in the table 1.

Simulation
Watchful reader noticed that table 1 does not include damping of pipe material and content of free air in water.To measure them is a big challenge, so they are subject to numerical optimization.
Numerical simulation with HYDRA software was done first.Damping and content of air were optimized to reach as good match of result with experiment as possible (especially, at the beginning of water hammer).Result in the figure 6 corresponds to b = 9.53 MPa s and M p = 2 • 10 −6 .First pressure surge is almost identical with experiment (sharp peak on the first pulse is higher).Frequency and amplitude are satisfying for the first three seconds, then simulation becomes less accurate.Probably, damping of pipe is not constant and some suitable function should be found or Voight-Kelvin model of solid should be replaced with some more complicated one.
Similarly, SimHydraulics simulation brought result in figure 7. Computed pressure does not match with experi-01037-p.3ment as well as previous simulation, but it is still sufficient.In this case, content of air and elasticity of the pipe were kept the same, just recalculated according to equations ( 4) and ( 7), but it was necessary to change value of damping to b = 12.925 MPa s.It, probably, caused very high peak at the beginning of water hammer (approximately 750 kPa).Frequency of pulsations is slightly higher than measured frequency, but optimization of damping, elasticity and content of air did not bring better result.
In both cases (HYDRA and SimHydraulics), pipe was split into 30 parts to ensure good accuracy of the computation.Greater number of segments did not bring noticeable difference.

Summary
Measurement of water hammer in gravitational pipeline was done.Static pressure in the pipe was about 100 kPa, so influence of free air in water has a big impact on the shape of pressure pulsations.Pulsations are also fast attenuated because of a strong damping in pipe material.It could be stressed that cavitation does not appear (the lowest pressure has value 17 kPa).
Computational simulation of water hammer was done with two pieces of different software: Matlab-Simulink-SimHydraulic product of Mathworks Company, and HY-DRA, product of Brno University of Technology.First one uses ode15s numerical method and the other one is based on Lax-Wendroff method.Both programs allow considering elasticity and damping of pipe and influence of free air in water, so they are supposed to solve presented case.
As value of pipe damping and quantity of free air could not be measured, they were optimized to find as good result of simulation as possible.
The best match with experiment was reached for volume fraction O p = 1.65 • 10 −3 at atmospheric pressure and damping b = 12.925 MPa s (according to SimHydraulics simulation) or 9.53 MPa s (according to HYDRA simulation).Difference could be caused by different numerical viscosity of both methods.

3 Fig. 1 .
Fig. 1.Dependence of density and bulk modulus of mixture wa-ter+air on the pressure, influence of poly-ethylene pipe elasticity is included.Mass fraction of air is 2•10 −6 and temperature is 15 • C

Fig. 2 .
Fig.2.Sound speed in the system according to fig.1

Fig. 3 .
Fig. 3. Response of pipe wall diameter on pressure change

Fig. 5 .
Fig. 5. Pressure in the vicinity of valve during water hammer

Fig. 7 .
Fig. 7. Comparison of SimHydraulics simulation (black line) with experiment (light line) both methods.Furthermore, Lax-Wendroff is more stable because of fixed time step.Elasticity and damping of pipe wall are defined directly by suitable values E and b.Again, Voight-Kelvin model of solid is used.HYDRA does not work with block scheme, but equations (2) and ( It has constant time step so solution takes longer time than ode15s, but only in case of steady state.When flow is nonstationary computational time is approximately same for 01037-p.2

Table 1 .
Important data of the experimental circuit Content of free air does not let the pressure drop under vaporizing value.Uncertainty of pressure sensor is ±0.25%.Sharp peak on top of the first pulse is caused by slow reaction of pipe wall when pressure increases.Strong damping of pipe material also quite fast attenuates pulsations.