Transient Simulations in Hydropower Stations Based on a Novel Turbine Boundary

Most accidents in hydropower stations happened during transient processes; thus, simulation of these processes is important for station design and safety operation. This study establishes a mathematical model of the transient process in hydropower stations and presents a new method to calculate the hydraulic turbine boundary based on an error function of the rotational speed. The mathematical derivation shows that the error function along the equal-opening characteristic curve is monotonic and has opposite signs at the two sides, which means that a unique solution exists to make the error function null. Thus, iteration of the transient simulation is unique and monotonous, which avoids iterative convergence or false solution and improves the solution efficiency compared with traditional methods. Simulation of an engineering case illustrates that the results obtained by the error function are reasonable. Then, the accuracy and feasibility of the mathematical model using the proposed solution are verified by comparison with model and field tests.


Introduction
Because hydropower stations play an important role in the peak regulation and valley filling of a power grid, the hydraulic turbine needs to frequently change its operating conditions and experiences many transient processes.The transient process in hydropower stations, including the interactions among hydraulics, mechanism, and electricity, is complicated.The closure of guide vanes and spherical valve induces a change in the flow inertia, which causes changes in the turbine rotational speed and hydraulic pressure in the piping system.When the working condition dramatically changes during transients, drastic changes in the waterhammer pressure and high rotational speed may lead to serious accidents that will endanger the safety of the hydraulic structure and turbine unit [1][2][3] and affect the power grid stability [4].Therefore, simulating the transient process of hydropower stations is necessary.The calculation accuracy is directly related to the design of the water diversion system, safe operation of the hydropower plant, and power quality.
The calculation methods of the water-hammer pressure are analytical [5,6], graphical [7,8], and computer-numerical [9].The analytical method is based on the chain equation of Allievi and the simplifications of the hydraulic turbine as a valve.This method is suitable for simple tube Pelton turbines.On the basis of the analytical method, the graphical method combines the water-hammer wave reflection and superposition with the turbine characteristic curves and guide-vane closing scheme.The traditional numerical simulation can be classified into two types: time-and frequency-domain methods.The time-domain method includes the method of characteristics (MOC) [8][9][10][11][12][13], finite difference (FD) [14,15], and finite volume (FV) [16,17].In these methods, the water delivery system is divided into a limited number of calculated cells.MOC is mostly applied owing to its advantages: high efficiency of computation, easy implementation for boundary conditions, easy parallelization with FD and FV for pipeline processing [18,19], and so on.By using MOC, various transition processes can be simulated.In the frequencydomain method [9], the basic equations are linearized to obtain the transfer function of the pipeline.Combined with the transfer functions of the governor and turbine, obtaining the dynamic responses of the hydropower system becomes easy.Irrespective of the time or frequency domain, the transient simulation model includes the pipe model and the turbine boundary.
The turbine boundary in the transient simulations can be divided into two categories [20].One is based on the geometric size of the turbine [21][22][23].The other is the characteristic curve based on the model tests, which is generally adopted in the transient simulations.The characteristic curves are usually transformed into curves with the unit speed as an independent variable and the unit discharge and unit torque as dependent variables.To enable transient simulation, the turbine characteristic curves are formulated into piecewise polygonal functions using an auxiliary grid [24][25][26].Then, the operating point consisting of unit parameters is expressed as By combining the two equations and other boundary equations, a one-element cubic equation of the square root of the head can be obtained, and the equation can be solved using the Newton-Simpson iterative method [27].The roots obtained by this method are correlated to the initial values and it can only obtain the root near the initial values.Meanwhile the first-order derivative of the function which is taken as the denominator should be neither too small nor null; otherwise, the iteration cannot be carried out.Generally, we take the root of last instant as the initial value.If the unit parameters obtained from the root of the iteration just lie in the line segment assumed previously, the root is considered as the correct solution and it will be taken as the initial values of the iteration in next instant.If the unit parameters are not in the assumed line segment, we should extend the search range and solve the cubic equation again until we get the right line segment or search all over the equalopening curve.If all of the line segments of the equal-opening curve cannot get a right root, we will reduce the accuracy of the iteration properly and search again.If all the methods proposed above cannot get a root, then we will select the root whose operating point is the closest to that of last instant as the right solution.Obviously, this calculation method suffers from two shortcomings.First, the root search direction is ambiguous; when the desired operating point is beyond the line segment, it is not clear whether the next search should target the last or next segment of the current line segment.Second, the roots obtained by the Newton-Simpson iterative method are correlated to the initial values; if the initial values are incorrect, the iteration results may be beyond the assumed line segment.Therefore, the equation may either become unsolvable or be incorrectly solved.This study establishes a mathematical model of a transient process by MOC and a novel turbine boundary, as shown in Figure 1 Two sets of characteristic curves are generally used to represent the turbine characteristics.Here, the traditional characteristic curves are grouped into a nonuniform spline surface [28], as shown in Figure 2. Thus, the turbine characteristics can be expressed as The measured data points (  1 ,   1 ,   1 ) of the turbine characteristics are first expressed as a grid -spline spacecurved surface, as shown in Figure 2(a).The space-curved surface is composed of three independent three-dimensional surfaces, as shown in Figure 2(b), or projections on a twodimensional plane, as shown in Figure 2(c).Parameter  represents the position of an operating point in the opening curve and is evaluated by uniform parameterization   = / ( = 0, 1, . . ., ).Parameter V stands for relative opening V  =   /  ( = 0, 1, . . ., ).The ranges of the two parameters are both in [0, 1], and each group of (, V) values corresponds to an operating point in the space-curved surface; that is, when  and V are determined, the three unit parameter values of the operating point become known.
According to the one-dimensional continuity and momentum equations of the pipe flow, the equations along the characteristic lines at the front and back of the turbine joint using MOC are as follows [14]:   ,   ,   , and   in (4) and ( 5) are all values of the previous instant and are known values in the current instant.  and   are determined by the flow and head of the previous instant, respectively, and   and   are both correlated to the wave velocity and pipe area, respectively, and are positive constants.
The equations for the unit parameters are expressed as The first derivative of the differential equation of the generator is expressed as Under an off-grid operation,   = 0 and   = 0. Thus, the generator equation can be simplified as The guide-vane opening is defined in advance; thus Parameter V represents the relative opening; therefore If the turbine is operating under a frequency or power control, ( 11) is eliminated, whereas the speed governor equation is added.
The unknowns in these equations include the heads   and   of the inlet and outlet of the rotating wheel, respectively, discharge   , rotating speed , torque   , unit speed   1 , unit discharge   1 , unit torque   1 , guide-vane opening , and parameters  and V.In general, eleven equations with eleven unknowns exist.

Calculation Steps of the Transient Process.
During the calculation of the load rejection in the transient process, the closure law of the guide vane is given, and the guide-vane opening is known at every instant.Thus, solving the abovementioned boundary conditions is simplified to seeking the operating points that satisfy the equations.With a known opening value, the equations can be solved as follows [27].
The flowchart is shown in Figure 3.When  1 > 0, the search interval of Ω = 0 is the whole equal-opening curve (i.e., from  = 0 to  = 1); when  1 < 0, the search interval is the part in which 1 is in the pump and reverse-pump regions.In general, this study needs to solve two problems: (1) searching for the operating point where Ω = 0 in a known equal-opening curve and (2) distinguishing between  1 and  2 when  1 < 0 and Assumption 1. Ω is a monotonic function in the search segment of the equal-opening curve with parameter .
Assumption 2. Ω has opposite signs at the two boundary points of this segment.
From Assumptions 1 and 2, a unique point must exist where Ω = 0.The search direction of the root can be determined according to the monotonicity, and the first problem in this study can be solved.For the second problem, that is, choosing between  1 or  2 , if Ω is a monotonic function, by inserting  1 or  2 into Ω, the value that can make the signs of Ω at the search boundaries opposite is required, and the other value should be discarded.In summary, to solve problems (1) and (2), we need to validate Assumptions 1 and 2.

Search Direction of 𝑋
We define Δ = √(Δ  1 ) 2 + (Δ  1 ) 2 + (Δ  1 ) 2 as the length of the segment, cos  = Δ  1 /Δ, cos  = Δ  1 /Δ, and cos  = Δ  1 /Δ as the cosines of the angles separately formed by segment    +1 and the three coordinate axes in space.Then Therefore, the sign of Ω/ depends on the sign of the sum of Δ  1 Ω/  1 , Δ  1 Ω/  1 , and Δ  1 Ω/  1 .When Ω/ > 0, Ω is a monotonically increasing function within    +1 ; when Ω/ < 0, Ω is monotonically decreasing.From the calculation method presented in Section 2.2 and ( 7)- (10), Ω can be expressed as The following parts will consider  =  1 as an example for analysis. =  2 can be analyzed by the same method.The partial derivatives of ( 17) are as follows: Evidently, Ω/  1 and Ω/  1 are constants that are greater than zero, and the sign of Ω/  1 varies under the following scenarios.(1) When the operating point is located in the pump region or pump braking regionwhere   1 < 0 and   1 > 0, Ω/  1 is positive.(2) When the operating point is located in the turbine braking region or reverse-pump region where   1 > 0 and   1 < 0, Ω/  1 is negative.(3) When the operating point is located in the turbine region and   1 < 2 4 1   1  (the estimated magnitude of   1 is not greater than 10 0 ), Ω/  1 is positive; otherwise, Ω/  1 is negative.The signs of Δ  1 , Δ  1 , and Δ  1 depend on the shape of the equal-opening curve.Because the characteristic curve of pump-turbine is the most complex, this study considers this curve as an example.In particular, the equal-opening curve of pump-turbine has the following characteristics [30] (as shown in Figure 4).Point A is the start point of the pump region.Point B is the inflection point of the pump operating region.Point C is the point where the unit discharge of the pump region is zero.Point D is the point where the unit rotating speed is zero.Point E is the point where   1 = 2 4  1   1  in the turbine region.Point F is the point where the unit torque is the largest in the operating region of the turbine.Point G is the point where the unit discharge is the largest in the operating region of the turbine.Point H is the point where the unit rotating speed is the highest in the operating region of the turbine.Point I is the runaway point.Point J is the point where the unit discharge is zero in the turbine region.Point K is the point where the unit rotating speed is the lowest in the braking region of the turbine.Point L is the end point of the reverse-pump region.

Comparison between
The magnitude of  4 1 is usually approximately 10 −4 , and the magnitude of  is usually in the range of 10 1 -10 2 .Hence, the magnitude of  4  1  is below 10 −2 .Therefore, along the entire equal-opening curve, the slopes approach infinity only in the BC segment of the pump region and at points H and K in the turbine region, as shown in Figure 5, where maybe  1 < According to ( 18) When   1 Ω/  1 and Δ  1 Ω/  1 are both greater than zero.Therefore, Ω/ > 0 along the characteristic curve, and Ω is a monotonically increasing function.
When  =  2 , this condition can be demonstrated in the same manner as  =  1 .Ω is a monotonically decreasing function in the pump region as well as in the reverse-pump region (where ).In summary, Ω is a monotonic function in the search segments; thus, Assumption 1 is verified.

Signs of Ω in the Search Boundary
) is on the left side of the   1 = 0 axis,   1,=0 <   1,= < 0, and   1,=0 <   1,= can be obtained.According to ( 7) and (14a), Ω =0 < 0 if  =0 >  = and  =0 <  = < 0. For the accordance, we take the root of last instant as the initial value.But when the right root is the further one, we will not obtain the right root based on the current initial value.Reducing the iterative accuracy or specifying one root according to last instant may result in bigger errors, and the unit parameters may deviate from the equal-opening curve and the transient parameters will jump up and down nonphysically.
The new method determines the unit parameters and the head at the same time.In the solving process, the equations of transient process and the operating point of the equalopening curve are mutually restrained.All roots obtained are consistent with the physical conditions in the calculation of the transition process: the operating points are in the equalopening curve; the rotational speeds calculated from speed equation and head equation are equal to each other.Therefore the present new solution method avoids the appearance of wrong roots.The results obtained are in line with the physical phenomenon without obvious vibrations.

Case Two (Model Validation).
A pumped-storage station model consisting of model units, piping systems, measuring system, and generator system was established in the laboratory.The schematic of the model is shown in Figure 9. Two different model pump-turbines, manufactured by Harbin Electric Corporation, were installed in the model station.The basic information of the pipe system and the pump-turbines is listed in Tables 3 and 4, respectively.
The basic information on the different sensors used in the experiments is presented in Zeng et al. [31].The positions of these sensors are shown in Figure 9. Various analog signals from these sensors were collected using an HBM Gen7i data acquisition system.The pump-turbine characteristic curve at the rated guide-vane opening was provided by Zeng et al. [31].
A load rejection experiment was conducted on the model station.The guide-vane opening of the pump-turbine was kept constant, and the experimental data are shown as the solid line in Figure 10.The simulation is also shown in  Figure 10 as a dotted line.The high consistency between the result of the experiment and the simulation confirms that the model has good accuracy.

Case Three (Field Validation).
A load rejection test of a pumped-storage power station in South China was conducted by ALSTOM.This section presents the comparison of the test results with the simulation results obtained through the mathematical model proposed in Sections 2 and 3.The power station has a similar layout with the case presented in Section 4.1.The basic information of the station is listed in Table 5.Two pressure sensors were installed on the unit.One was located downstream of the ball valve 3 m away from its center line.The other was located at the ell of the draft tube.The altitude of this sensor was 133.12 m.By comparing the results of field test and simulation, we can clearly see and thus could arrive at the conclusion that the simulation has good accuracy and can reflect the hydraulic transients of the load rejection with acceptable errors.These errors may be caused by the inaccurate characteristic curve, imprecise parameters of the piping system, and errors in the installation location of the sensors and cross-sectional area of the tubes.In addition, the one-dimensional characteristic method cannot simulate the pressure fluctuation, as clearly shown in Figures 11(b) and 11(c) at the peak of the spiral case pressure and nadir of the draft tube pressure, leading to inconsistency between the field test and the simulation.

Conclusions
This study has presented a mathematical model of the hydraulic transient process and proposed a new method to solve the turbine boundary.The new solution was based on error function Ω of the rotational speed.According to the change in the values and signs of the error function along the equal-opening characteristic curve, we can discriminate and search the root of the transient equations; when Ω = 0, the iteration in the numerical solution was completed.This study has shown that Ω is a monotonic function in the search segment along the equal-opening characteristic curve and that its signs in the search boundary are opposite to each other.The sign of Ω may be used to determine the search direction of the operating point in the equal-opening characteristic curve in the next time step.Within the search range, only one point exists that could fulfill the condition Ω = 0; hence, the equation set has only one unique solution.

Mathematical Problems in Engineering
The simulation of the engineering case illustrates that the operating trajectories obtained by the method proposed in this paper are more reasonable and credible, and the pressure at the end of the spiral case and at the entrance of the draft tube obtained by this method does not sharply jump.In contrast with the model and the field test results, the simulation results can accurately reflect the change process of the transient parameters in the load rejection.Therefore, this method is clearly superior to the traditional solution and can improve the simulation accuracy of the transient process of a hydropower station.

Figure 3 :
Figure 3: Flowchart of the solution process of the turbine boundary equations based on the space surface.

3. 1 .
Directional Derivative of Ω.According to the directional derivative theorem, the directional derivative of Ω in line segment    +1 shown in Figure2can be expressed as[29]

Figure 4 :
Figure 4: Key points of the characteristic curves.

Figure 8 :
Figure 8: Comparison of the results of the traditional and new solutions.

Figure 10 :
Figure 10: Comparison of the simulated parameters and experimental data.

Figure 11 :
Figure 11: Comparison of the simulation and field test.
22 41 , both  1 and  2 are positive real numbers, and we need to select the right one.Under other conditions, when both  1 and  2 are nonpositive real numbers, Step (1) is repeated.

Table 1 :
Signs of the different terms at different segments.

Table 3 :
Basic unit information of the pipe system.

Table 4 :
Basic unit information of the pump-turbines.

Table 5 :
Basic unit information of the pump-turbines.