The Possibilities of Dynamic Phenomena Calculation in the Penstock of Hydropower Plant via Different Methods

. When designing a hydropower plant with a long pressure penstock, special attention must be paid to the pressure conditions in the penstock, especially during emergency flow shutdown. In this paper, the main principles of the calculation of static pressure rise in the pressure pipeline during rapid flow shutdown are presented. The results obtained by analytical methods (calculation of direct and indirect water hammer) and numerical methods (mathematical modelling in Matlab software environment) are compared on specific examples. The influence of some input parameters of the numerical model on the results of the numerical calculation is analysed separately.

Among these phenomena we can also include the pressure increase in the penstock caused by the rapid shut-off of the flow at the end of the penstock. The reason for the rapid flow closing may be e.g. an emergency condition on the turbine or a sudden de-phasing of the generator accompanied by a subsequent increase in the turbine speed up to the value of runaway speed. In these circumstances, the water supply to the turbine must be stopped urgently (by closing the distributor or the shut-off valve in front of the turbine). In this case, however, there is (theoretically) a risk of a water hammer occurring at very rapid closure, which would have destructive effects on the penstock or the other equipment. It is therefore essential to set the emergency shutdown time of the equipment appropriately so that the maximum pressure rise in the penstock or upstream of the turbine does not pose a safety risk. Appropriate closing time setting requires calculating the pressure rise for several alternative closing times and selecting a value that will not cause a dangerous pressure rise and at the same time will not cause an unnecessarily long duration of the turbine runaway speed.
The above problem can be solved by several different methods. One possibility is to calculate the pressure increase with the application of classical theory and analytical (explicit) equations applying specific idealizations of the real system. The other option is to use mathematical modelling and numerical mathematics methods that are now common in commercial software such as Matlab & Simulink.
The aim of this paper is to compare the calculation of the pressure rise during penstock closure using analytical methods (idealized explicit equations) and numerical methods. At the same time, the influence of the time step size on the results obtained using numerical calculation methods is analyzed. The comparison is made for a simple piping system (constant diameter along the length, homogeneous material and constant wall thickness) which can form the penstock of a hydropower plant.
The phenomenon of water hammer became an important research topic in the 19th century for the field of water turbines and pumps. The experimental work of Zhukovsky (1900) defined the foundations of the current theory of water hammer. His research enabled the formulation of the fundamental equation in the theory of water hammer, which relates the change in pressure to the change in velocity at instantaneous valve closure [5].
Allievi (1902) developed a mathematical approach that analytically solved the piping system and found pressure values for different time steps and locations in a one-dimensional model of a thin-walled pipe for any linear flow variation. In 1913, Allievi presented a general theory of water hammer based on his previous mathematical approach. The applicability of the equations to various model problems was achieved by introducing a characteristic parameter , which describes a certain element of the pipeline [3]. Jaeger (1933) developed the general theory of Allievi for specific cases. He developed analytical solution methods for systems with more than one pipe, thick-walled pipes, concrete-covered pipes, and systems with equalization chambers. The analytical methods based on Allievi's theory provide a general understanding of the water hammer phenomenon, but the formulas derived from it are limited to only a few model cases. The methods neglect friction and the second-order terms contained in the system of hyperbolic partial differential equations. They also do not allow the solution of problems with multiple characteristic parameters . They provide a solution only for the pressure field and treat the flow drop during valve closure as linear [3].
Due to the existence of certain limitations of the analytical solution methods, so-called graphical solution methods of the hydraulic order have been introduced. Schnyder's (1936) and Bergeron's (1935) graphical methods gave results for the pressure and velocity fields for each time step at any location in a one-dimensional model of a thin-walled pipe for any type of shut-off valve. Parmakian (1963) presented a modified Schnyder-Bergeron graphical method that is applicable to systems with multiple characteristic parameters , with consideration of hydraulic losses and with composite piping. The modified graphical method provides an accurate and fast solution for simple piping systems. However, for complex piping systems containing many hydraulic elements and branches, the graphical method did not prove successful and was replaced in a short time. Numerical methods [3], [6] solved the above problem.

Numerical solution methods
The development of computer technology has made it possible to solve hyperbolic partial differential equations by numerical methods. Gray (1953) used the Eulerian method of characteristics (MOC). The method allows one to determine the order of approximation while neglecting the effects of friction, kinetic energy, and pressure wave shape change. Lister (1960) compared two methods of MOC discretization: the method of specified time intervals also called the fixed grid method (Fig. 2) and the so-called characteristic grid method. He found that the fixed grid method is much simpler in terms of computation and also in terms of constructing the computational grid. Because of its accuracy, simplicity, numerical efficiency and programming inexpensiveness, the MOC discretized using the fixed grid scheme became the most widely used method for solving the water hammer in its time. Perkins et al. (1964) established a numerical stability condition: / ≤ 1/ For the characteristic fixed-grid scheme, which is also referred to as the Courant stability condition. Lai (1961) in his work included a nonlinear hydraulic loss term in the MOC and implemented it in a computer code [6]. The Lagrangian method, also known as the Wave Plan Method, was introduced by Wood et al. (Wood later compared WPM with MOC and found that for tasks such as large pipelined networks, the wave plan method requires less CPU time than MOC. Chaudhry (1985) proved that the second order finite difference method (FDM) gives better results than the first order MOC. However, he found that this method increases computation time and data storage size requirements. The first-order finite volume method (FVM) of Gudonov type was first used to solve water hammer problems by Guinot (2002). Zhao (2004) developed an explicit second-order Gudonov type method and compared its results with MOC with spatial line interpolation. He found that the second-order Gudonov type explicit finite volume method achieves a given level of accuracy with significantly less CPU time than MOC with spatial line interpolation [8]. In order to predict the pressure change caused by valve closure in a gravity pipe with air present, a numerical model based on the Lax-Wendroff method was developed. Compared to traditional methods, this study provides a more advanced approach by taking into account the influence of the flow velocity in the pipe on the propagation of the pressure wave and on the simulation of transient processes. In practice, there are many different water pipeline systems operated with the presence of a certain air content, therefore the Lax-Wendroff Method (LWM) is significant for its relatively accurate prediction of water hammer processes, especially under these conditions [7]. Fig. 3 compares the processed results from the simulation of water hammer by HYDRA software which uses the explicit Lax-Wendroff

Analytical calculation of the pressure rise in the penstock
Analytical calculation methods apply equations to solve for the impact pressure in the penstock. A more detailed derivation and explanation of the calculation procedure is given, e.g., in [2]. For a penstock with a given length, constant diameter and wall thickness, equation (1) holds for the reflection time of the shock wave. Where: -the bulk modulus of elasticity of the liquid, -liquid density, -internal diameter of the pipe, -modulus of elasticity of the pipe wall material, -pipe wall thickness.
In the event of sudden closure of the pipe (closure time is less than the reflection time of the shock wave), a direct water hammer occurs. Then (3) applies.
The increase in static pressure after sudden closure of the pipe (in direct water hammer) can be determined according to equation (4) or the increase in pressure can be expressed in metres of water column (5). Where: 0 -steady flow velocity of the liquid in the pipe (before closing or before the start of closing).
In the case of a gradual closure of the pipe, where the closure time is longer than the reflection time of the shock wave (6), an indirect water hammer occurs. The increase in static pressure, expressed in metres of water column, can be determined by equation (7), which applies to a 'loose' pipe. The above calculation procedure can be applied to calculate the pressure rise in the model penstock according to the scheme in Fig. 4. The pre-closure flow rate, gradient, penstock length, pipe diameter, wall thickness and other parameters required to calculate the pressure rise are given in Table 1. For the case of a direct water hammer, an almost instantaneous pipe closure (0.05 s) is considered, which is only a hypothetical value that is in fact technically difficult to achieve, and for an indirect water hammer, a closure time of 12.3 s is considered. The calculated pressure surges for the closure times considered are also given in Table 1.

Numerical simulation of water flow
The analytical models allow the fast calculation of the parameters needed to determine the main design variables of a given work. However, it is more appropriate to use numerical methods and simulation programs to calculate all the variables, or to calculate them with consideration of time changes. A suitable tool is, for example, the commercial software Matlab or its open-source version SciLab. These tools allow 0D/1D and 2D analysis of physical problems and tracking of selected physical quantities as a function of changing boundary conditions and parameters.
For the purpose of the given simulation, Matlab software with Simulink superstructure was used, which simplifies the calculation by representing the steps and dependencies of the computational simulation by function blocks.   5 shows the blocks used to set the boundary conditions. The "computation setup" block defines the basic parameters of the computation, i.e. mainly the solver type, the number of iterations and the time step of the computation. The "fluid properties" block defines the type of fluid and the properties to be considered by the calculation. The ideal reservoir block defines the reservoir type and the surface pressure. The calculation blocks themselves represent real parts of the technical equipment in question, i.e. the penstock pipe and the shutoff fitting, the valve. For simulation purposes, a pressure sensor and a flowmeter are still included in the network to record the actual state of these quantities.   Fig. 6. From Fig. 6 it can be seen at which values of time step approximately the same simulation results are obtained.
The actual results of the numerical calculation of pressure rise performed for the closing time of 12.31 s and comparison with the result of the analytical calculation is shown here in Table 2.

Conclusion
The paper presents the results of analytical and numerical calculations of the pressure rise in the penstock of a hydropower plant caused by emergency shutdown the flow. The results of the analysis of the influence of the time step on the obtained results of the numerical calculation are presented. The comparison of the analytical and numerical calculation results showed that the result of the numerical calculation of the pressure rise in the penstock is more than 10 % different. For a more detailed evaluation of the individual calculation methods accuracy, it would be necessary to compare the obtained results with the results of operational measurements on the penstock. It should also be noted that there is still a relatively large potential for further improvement of the mathematical model applied in the numerical calculations.
Authors would like to thank to project: Regenerácia použitých batérií z elektromobilov / Regeneration of used batteries from electric vehicles ITMS2014+: 313012BUN5