Computational compliance criteria in water hammer modelling

Among many numerical methods (finite: difference, element, volume etc.) used to solve the system of partial differential equations describing unsteady pipe flow, the method of characteristics (MOC) is most appreciated. With its help, it is possible to examine the effect of numerical discretisation carried over the pipe length. It was noticed, based on the tests performed in this study, that convergence of the calculation results occurred on a rectangular grid with the division of each pipe of the analysed system into at least 10 elements. Therefore, it is advisable to introduce computational compliance criteria (CCC), which will be responsible for optimal discretisation of the examined system. The results of this study, based on the assumption of various values of the Courant-Friedrichs-Levy (CFL) number, indicate also that the CFL number should be equal to one for optimum computational results. Application of the CCC criterion to own written and commercial computer programmes based on the method of characteristics will guarantee fast simulations and the necessary computational coherence.


Introduction
When water flow is started or stopped faster than the system can respond acoustically, a pressure surge commonly referred to as a water hammer event will occur.The same principle is also applicable when the flowing medium is another fluid (oil, steam etc., in the case of steam, the pressure surge is known as a steam hammer event).In the mentioned transient state, water hammer waves -also called pressure waves -are transmitted through the system at the speed of sound of the water-pipe combination.In typical steel pipes, the pressure wave speed is greater than 1200 m/s.This pressure waves can cause major problems, from noise and vibration to pipe collapse and total system failure.This is the reason why the protection of water systems or nuclear power plants cooling systems against water hammer is so important.Designers of such a pressurised systems normally use expensive additional equipment to protect them against water hammer damage: hydropneumatic or surge tanks, special valves, air chambers, rapture disks [1,2].But the simplest, effective and nearly cost free protection is to control the valve-closing time and the pump speed.For optimal control procedures at the stage of designing such systems, a numerical simulation analysis is needed [3,4].It is assumed that the flow of this type is being described by the hyperbolic system of partial differential equations (sequentially by the continuity equation and the motion equation) [1][2][3][4][5][6]:   In this connection, the stability criterion is not responsible for the accuracy of modelling but only for its correctness in the mathematical sense.Following the analysis of the method of characteristics (which is used to numerically solve the partial equations that describe the water hammer), as well as the results being obtained themselves, it is reasonable therefore to introduce another criterion, the objective of which would be to guarantee compliance between simulated results and those being observed experimentally.

Purpose and scope of the study
Unfortunately, the studies that used the method of characteristics to model the unsteady pipe flows have not indicated at what minimum division the results being obtained will be characterised by an acceptable error.Simpson and Bergant [8] modelled the pipe flows with cavitation assuming different computational elements (Table 1).They observed, when applying the division into 32 elements and more, that there are significant deviations and a significant variation in the results being obtained.This fact is explained by them that the greater the number of divisions into grid reach lengths, the more cavitation areas is modelled and, at the same time, the greater is the likelihood of a superposition of waves.They suggested that the number of reaches (grid points) over the pipe length (N) is to be chosen so that the ratio of the maximum cavitation volume being modelled to the volume of fluid contained in a single numerical element is to reach the value below 10%.Interestingly, they did not notice any modelling problems when applying division into 4 or 8 elements.In 2001, in another study -this time on the pipe flows without cavitation-related column separation, Bergant et al. [9] applied the division only into N = 8, 16 and 32 elements and showed that different number N on the diamond (staggered) grid, using the friction model according to Brunone, does not affect the quality of simulated pressure runs (Fig. 1).-the method of characteristics is less sensitive to the change of N factor than other methods (for example, the finite-difference method -the studies by Hadj-Tajeb [13,14]).On the basis of the literature being reviewed above, the objective and scope of the study could be defined.The purpose of tests performed and described in this study was to determine the optimum value (number) of the numerical division of hydraulic pipe (of the sample length L) into computational elements (marked in this study as N).This division should be as small as possible so that the computer calculations that use the method of characteristics are to be made as fast as possible.
Owing to the fact that the rectangular grid of the characteristics is commonly used and the effect of N on them was analysed according to Commonly, a sharp toothed increase in pressure is being observed on the peak of the first pressure amplitude on the rectangular grid.This drawback can be minimised, however, by showing on the graph illustrating the variation of pressure in a given pipe cross-section being analysed every second simulated result.Then, a smoothed simulation waveform will be obtained, very similar to the results being obtained on the diamond grid.
The scope of the study will refer to the comparative analysis of a number of curses being modelled (with and without cavitation; laminar, transient and turbulent: 1100<Rei<40350) with the experimental ones being obtained in a typical hydraulic system.

Effect of number of reaches on simulation results
In the tests presented below, the hydraulic pipe being examined was divided into N computational elements, starting from N = 2, every 1, to N = 64.Carrying out computer simulation in such a wide range guaranteed an accurate picture of the variation of two basic quantitative parameters: Et and Ep, being determined for each simulated pressure run.Parameters Ep (maximum pressure compliance) and Et (time compliance) characterise compatibility of the pressure runs being modelled [19].They are calculated using the experimental and numerical results (Fig. 6) from the following formulas: where: k -number of analysed pressure amplitudes; pis and pie -maximum pressure on the i-th analysed amplitude from the simulated and experimental graph, respectively; tis and tie -time of occurrence of the maximum pressure on the ith analysed amplitude from the simulated and experimental graph, respectively.Fig. 6.Analysed pressure and times values in the quantitative analysis.

In his study, when determining parameter Et, the time compliance on the first amplitudes was omitted because the chosen experimental results used for comparisons were characterised by a noticeable effect of the vibrations of the valve being closed on the initial phase of pressure curve (fluid structure interaction (FSI) effect). Experimental tests necessary for comparative analysis have been performed by Adamkowski and Lewandowski at the Institute of Fluid-Flow Machinery of the Polish Academy of Sciences in Gdańsk. In their paper of 2006
[11], these authors showed courses, in which cavitationrelated column separation did not occur.From that study, nine courses were selected for the comparisons being made in this study (Table 2).However, in the study of 2012 [12], the same authors showed different results obtained on a slightly modified test bench, in which cavitation-related column separation occurred.In the comparative tests presented below, six courses obtained then were used (Table 3).The details necessary for numerical modelling are compiled in the presented tables (Table 2 and 3).The values of pressure wave propagation speed c, which were used in computations, collected in the Tables 2 and 3 were not calculated from theoretical equations but determined by the analysis of experimental pressure runs.

Fig. 7. Results for Ep (cases without cavitation).
The Table 2 shows that the average value of pressure near the reservoir pT was 1.263•10 6 [Pa] and that the differences observed in the speed of pressure wave propagation (1304-1308 m/s) were slight (from the average value of 1305.44 m/s, the error was only ± 0.2%).For the tests with cavitation-related column separation (Table 3), it can be concluded that the differences being observed experimentally in the speed of pressure wave propagation were significant in relation to those without cavitation (from the average value of 1246.667m/s, the error was ± 1.36%).In the tests being performed, dispersion of the values of pressure wave propagation speed c was as much as 30 m/s, suggesting that there are bubbles of air in the flowing fluid in the chosen pressure runs, in which the propagation speed was the lowest because this is mainly them that decrease the speed of pressure wave propagation.with cavitation their number varied depending on the boundary and initial conditions (for k values see Table 3).On the other hand, the next two graphs show the variation of the same parameters but for the courses in which column separation occurred, i.e. there were cavitation areas appearing (Fig. 9 and Fig. 10).From the results of quantitative tests shown in a graphical form (Fig. 7 to Fig. 10), the following conclusions can be drawn: • in the pressure runs without cavitation, the acceptable compliance of simulations with the experimental results was obtained for the division into N = 10 (Fig. 7 and Fig. 8).Analysis of the variation of parameter Ep (Fig. 7) shows that this error will be smaller for N > 60.On the other hand, analysis of the variation of parameter Et (Fig. 8) showed that with the increasing N in five runs (v0 = 0.162; 0.34; 0.467; 0.559; 0.631) there was a slight but steady increase in the error.However, it can be considered that the Et results being obtained for N = 10 were characterised by acceptable values; • analysis of the results obtained for courses with cavitation (Fig. 9 and Fig. 10) showed that in these runs the results can be accepted if division of the pipe is assumed at 10 elements (N = 10) although smaller dispersion of the results in these analysed runs is visible for N > 20; • pressure runs with cavitation in the transitional range (v0 = 0.47; 0.82) are characterised by unacceptable modelling errors Ep within the ranges of 8% to 16% and a high sensitivity to changes in parameter N (Fig. 9).This may suggest the need of further work on improving the applied cavitation model for the range of transitional pipe flows.

Effect of the Courant-Friedrichs-Lewy (CFL) number
The Courant-Friedrichs-Lewy (CFL) number expresses the ratio of the distance traveled by a disturbance in one time step to the length of a computational distance step.As pointed in many textbooks, the CFL number in the method of characteristics should be less than or equal to unity [1, 2, 5, 6] so as to ensure that the solution remains within the computational domain.However, the exemplary numerical simulations performed in this study, consisting in the application of the CFL number from the following range 0.5<CFL≤1 (the following values were adopted for the tests: CFL = 1, 0.95, 0.9, 0.85, 0.8, 0.7, 0.6 and 0.5) for the same simple system, showed that the application of lower value of the CFL number results only in deterioration of the modelling

Conclusions
A number of simulations performed in this study helped to determine a numerical criterion of the compliance of simulated results.The computation compliance criterion CCC enforces at the beginning of numerical calculations the necessary division of pipe into computational elements.According to the test results being observed in the method of characteristics for rectangular grids, the following condition must be satisfied: An additional condition, which was demonstrated in section 3, is to maintain the stability criterion in the form of the Courant-Friedrichs-Lewy number equal to unity (CFL = 1).In order to confirm the condition defined above in a wider range, further tests for pipes with other lengths, as well as other internal diameters, are necessary.In addition, the performed tests clearly show that the unsteady pressure runs without cavitation are simulated with a greater accuracy; this may suggest that there is a need for modification of the cavitation model according to Adamkowski-Lewandowski so that the error for the Reynolds number, especially from the range of transitional flow (multiple transition from the laminar flow in the turbulent one during the pressure pulsation dumping is a phenomenon that accompanies the water hammer event) is also below the acceptable range of 5%.

Fig. 1 .
Fig. 1.Bergant et al.'s results for different N [9].Bughazem and Anderson [10] pointed to a significant effect of N on simulation results based on a series of tests performed on the rectangular grid of the characteristics (they used the friction model according to Brunone).Adamkowski and Lewandowski, when examining the effect of friction on simulation results of the pressure runs without cavitation [11], observed no effect of the grid density on the simulation compliance and the stability of numerical computations when applying the division into N = 12, 24 and 36.In their next comparative study on the pipe flows with cavitation [12], the same authors showed that the division of N, however, affects the time of modelling the transition from the course with cavitation into the course without it (Fig. 2).

Fig. 2 .Fig. 3 .
Fig. 2. Envelope of the analysed pressure curve obtained for different N [12].The significant effect of N on simulation results (Fig.3, Fig.4) has been also observed by Hadj-Taieb et al.[13,14], but these authors did not use the method of characteristics and solved the basic equations with the finite-difference method (Lax-Wendroff method).However, while analysing the results obtained by Karadžić et al.[15], it is clear that also for the relatively dense division N = 54, 108 and 216 significant difference between the simulated results can be seen (Fig.5).In other studies, their authors often forget to specify which numerical discretisation was used to perform comparative tests[17][18][19][20] or they only show the applied

Fig. 5 .
Fig. 5. Variation of results for large N [15].Summing up, the following three important conclusions can be drawn from the studies being analysed: -no significant effect of N on simulation results is being noticed in the pressure runs without cavitation [9, 11]; -the effect of N seems to be much larger in the pressure runs with cavitation, although this does not refer to the modelling of the first amplitudes being characterised by maximum pressures but rather to the time of transition of the courses with cavitation-related column separation into the courses without cavitation [8, 12, 15]; Owing to the fact that this study examined the effect of N on the pressure runs with cavitation, it was necessary to implement additional conditions for determination of the volumes of emerging cavitation areas in the computer programme being used to solve the basic equations describing the unsteady flow of that type.As a model of cavitation-related damage of the flow continuity, the model according to Adamkowski-Lewandowski was used, which has been described in detail in the paper of 2009[27].

: ∆t -time step [s], ∆x -single reach length (horizontal distance between neighbouring grid nodes) [m]. The exemplary simulations performed in this study indicate that this criterion should be rather written in the form of equality. As for the compliance criterion, it is difficult to find in the current literature on this subject any condition which would have to be satisfied to maintain necessary simulation accuracy. As shown by the results of the comparative tests performed and
E3S Web of Conferences 19, 03021 (2017)