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.

1. 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 water-hammer pressure and high rotational speed may lead to serious accidents that will endanger the safety of the hydraulic structure and turbine unit [13] 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) [813], 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 frequency-domain 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 [2123]. 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 [2426]. 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 equal-opening 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. In terms of the novel turbine boundary, an error function on the rotational speed is constructed, and a new solution method based on the monotonicity of this error function is presented. The novel solution method is theoretically analyzed, and the superiority of the new solution is validated by transient simulations, model test, and field test.

2. Mathematical Model and Calculation Method

2.1. Mathematical Model of a Novel Turbine Boundary

The mathematical model of the transient process comprises the pipe, turbine boundary, generator, and speed governor models.

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 of the turbine characteristics are first expressed as a grid -spline space-curved 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 two-dimensional 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 . Parameter stands for relative opening ). The ranges of the two parameters are both in , and each group of values corresponds to an operating point in the space-curved surface; that is, when and 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, and . Thus, the generator equation can be simplified as

The guide-vane opening is defined in advance; thus

Parameter 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 , unit discharge , unit torque , guide-vane opening , and parameters and . In general, eleven equations with eleven unknowns exist.

2.2. 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].

(1) We assume a value for parameter and obtain the three unit parameters from (3). Inserting , , , and into (4)–(6) and simplifying the result yield

Then, the two roots are obtained as follows: represents the square root of the head and is always a positive real number. According to (14a) and (14b), when . When and , both and are positive real numbers, and we need to select the right one. Under other conditions, when both and are nonpositive real numbers, Step is repeated.

(2) We insert into (7) and (8) and obtain rotating speed and torque . We insert into (10) and obtain rotating speed .

(3) We define . When , parameter is reasonably assumed. and the three unit parameter values of the operating point are desired. Hence, we proceed to Step (4); otherwise, Steps (1)–(3) are repeated.

(4) We insert into (4)–(6) to obtain the values of , , and .

The flowchart is shown in Figure 3. When , the search interval of is the whole equal-opening curve (i.e., from to ); when , the search interval is the part in which is in the pump and reverse-pump regions. In general, this study needs to solve two problems: searching for the operating point where in a known equal-opening curve and distinguishing between and when 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 . 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 or , if Ω is a monotonic function, by inserting or 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 and , we need to validate Assumptions 1 and 2.

3. Search Direction of

3.1. Directional Derivative of

According to the directional derivative theorem, the directional derivative of Ω in line segment shown in Figure 2 can be expressed as [29]

We define as the length of the segment, , , and as the cosines of the angles separately formed by segment and the three coordinate axes in space. Then

Therefore, the sign of depends on the sign of the sum of , , and . When , Ω is a monotonically increasing function within ; when , Ω is monotonically decreasing. From the calculation method presented in Section 2.2 and (7)–(10), Ω can be expressed as

The following parts will consider as an example for analysis. can be analyzed by the same method. The partial derivatives of (17) are as follows:

Evidently, and are constants that are greater than zero, and the sign of varies under the following scenarios. When the operating point is located in the pump region or pump braking regionwhere and , is positive. When the operating point is located in the turbine braking region or reverse-pump region where and , is negative. When the operating point is located in the turbine region and (the estimated magnitude of is not greater than 100), is positive; otherwise, is negative.

The signs of , , and 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 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.

Table 1 lists the changes in the signs of , , and along the equal-opening curve.

3.2. Determination of the Monotonicity of
3.2.1. Comparison between and

Assuming as the ratio between and , we letAccording to (18)The magnitude of is usually approximately 10−4, and the magnitude of is usually in the range of 101–102. Hence, the magnitude of 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 , ; however, these areas are also the segments where . According to the signs of these three terms in Table 1, we know that . For the other intervals of the equal-opening characteristic curve, .

Consequently, except for the BC segment and the neighborhoods of points H and K, the sign of depends mainly on the signs of and , whereas can be ignored.

3.2.2. Comparison between and

Assuming as the ratio between and , we letAccording to (18)When , .

Figure 6 shows how the signs of and change along the equal-opening curve. is determined by the shapes of the characteristic curve. It reaches a low value in the BC segment and at points H and K in the entire opening curve while reaching a peak at point D. reaches its maximum value when and reaches its minimum value near the axis. In the pump braking and turbine operating regions where smoothly changes (the range of CG or CH, including EG), . At the ends of the pump and reverse-pump regions, is also possible; however, the two positions are usually beyond the calculation range of the transient process. At the other positions, exists.

When , opposite signs for and occur at two positions. The first position is the curve segment EG where , , and . Then, . The second position is the curve segment HK where , , and . Then, . At other positions, and are both greater than zero. Therefore, along the characteristic curve, and Ω is a monotonically increasing function.

When , this condition can be demonstrated in the same manner as . Ω 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.

3.3. Signs of in the Search Boundary
3.3.1. When

If the current operating point is on the left side of the axis, , and can be obtained. According to (7) and (14a), if and . For the end point in the reverse-pump region, and if . As a result, when , Ω is a monotonic function within the range of and . Therefore, only one point exists where along the entire characteristic curve. The same conclusion can also be drawn when the operating point is on the right side of the axis.

3.3.2. When

Similarly, when , and . When , and .

If , for the boundary point where in the pump region, we can obtain the result by inserting into (14a) and (14b). Then, inserting into (17) yields the following result:At the position where is relatively larger, the second and fourth terms in (23) are smaller than the other two terms; thus, they can be ignored. By comparing the first term with the third term, when ; otherwise, . Similarly, at the boundary point where in the reverse-pump region, when ; otherwise, .

In summary, when , (or ) in the search boundary of the pump region (or reverse-pump region) only if . When , (or ) in the search boundary of the pump region (or reverse-pump region) only if .

Based on the systematic discussions of and , two conclusions can be drawn. The signs of Ω at two boundary points in the search segment are opposite; thus, Assumption 2 is verified. When , ; meanwhile, when , . Thus, the second problem is solved.

4. Engineering Cases

4.1. Case One (Theoretical Validation)

To validate the proposed solution, simulations of the load rejection set in a pumped-storage power station were carried out and presented in this section. The details of the pumped-storage power station, along with simulation results of both the traditional method and proposed solution, are presented as follows.

The simulated pumped-storage power station has four units (300 MW each and 1200 MW in total) supplied by a single tunnel, and two surge chambers were installed on both upstream and downstream. Figure 7 shows the schematic layout of the station. Table 2 lists the basic information of the single machines and the parameters of the calculated case.

During the simulation, all turbines simultaneously rejected their loads, except for one that experienced a guide-vane mechanical failure and its opening remained unchanged. For the other three turbines with no malfunction, their guide vanes held their position for 10 s and then completely shut in 15 s. Figures 8(a)8(d) show the difference between the transient parameters using the new solution proposed in this paper and the traditional solution.

Until 20 s, the operating trajectories and variations in the flow and pressure at the inlet and outlet of the rotating wheel obtained by traditional method are more logical and credible with no obvious fluctuation. After 20 s, the operational trajectories become thicker at the reverse-S-shaped region of the turbine characteristic curve and even deviate from the scheduled equal-opening curve. Obvious vibrations occur in both the flow and pressure values, indicating that the traditional method is unable to determine the correct operating points. In contrast, no obvious vibration is observed in the calculation results obtained by the proposed method in this paper, and the calculation results until 20 s coincide with those obtained by the traditional solution. Therefore, this new method is clearly superior to the traditional solution and can improve the calculation accuracy of the transient process.

The traditional methods is to solve the head first and then seek the operating point. There may be more than one solutions to the one-element cubic equation and the roots obtained by the Newton–Simpson iterative method are correlated to the initial values. Because lack of sufficient 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 equal-opening 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 equal-opening 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.

4.2. 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.

4.3. 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.

5. 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 , 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 ; hence, the equation set has only one unique solution.

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.


:Interpolated coefficients of characteristics
: Interpolated coefficients of characteristics
:Diameter of runner inlet [m]
:Self-regulation coefficient of power network
:Flywheel moment [kg·m2]
:Pressure head [m]
:Rotating inertia (= ) [kg·m·s2]
:Resistance moment of generator [kg·m]
:Driving force moment of hydraulic turbine [kg·m]
:Unit torque [kg·m]
:Specific speed (m·kW)
:Rotational speed [r/min]
:Unit speed [r/min]
: Rated power [MW]
:Cross-section discharge [m3/s]
:Unit discharge [m3/s]
:Time [s]
:Position of an operating point in opening curve
:Relative opening degree
:Square root of head []
:Absolute guide-vane opening [mm] or [°]
:Angular velocity (= ) [rad/s].
:Value of last time step
:Generator parameter
:Maximum value
:Rated value
:Head of the runner outlet
:Hydraulic turbine parameter
:Head of the runner inlet.

Competing Interests

The authors declare that there are no competing interests regarding the publication of this paper.

Authors’ Contributions

Yanna Liu developed the theory, produced the results, and wrote the paper; Jiandong Yang conceived the project, provided measurement data of hydropower plants, and supplied guidance as a supervisor; Jiebin Yang performed part of programming work and conducted part of case studies and discussions; Chao Wang performed part of programming work and gave technical support; Wei Zeng performed the model test and polished the paper.


This work was supported under the Key Program of the National Natural Science Foundation of China (Grant no. 51039005).