# Transient Simulations in Hydropower Stations Based on a Novel Turbine Boundary.

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 [1-3] and affect the power grid stability . 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 . 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-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 frequency-domain method , 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 . One is based on the geometric size of the turbine [21-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-26]. Then, the operating point consisting of unit parameters is expressed as

[Q'.sub.1] = [f.sub.1] ([n'.sub.1], y) = [A.sub.1] + [A.sub.2][n'.sub.1] (1)

[M'.sub.1] = [f.sub.2] ([n'.sub.1], y) = [B.sub.1] + [B.sub.2][n'.sub.1]. (2)

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 . 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 maybe 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 B-spline surface , as shown in Figure 2. Thus, the turbine characteristics can be expressed as

[mathematical expression not reproducible]. (3)

The measured data points ([n'.sub.1], [Q'.sub.1], [M'.sub.1]) of the turbine characteristics are first expressed as a grid B-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 u represents the position of an operating point in the opening curve and is evaluated by uniform parameterization [u.sub.i] = i/m (i = 0, 1, ..., m). Parameter v stands for relative opening [v.sub.j] = [y.sub.j]/[y.sub.n] (j = 0, 1, ..., n). The ranges of the two parameters are both in [0,1], and each group of (u, v) values corresponds to an operating point in the space-curved surface; that is, when u 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 :

[Q.sub.P] = [Q.sub.CP] - [C.sub.QP] x [H.sub.P], (4)

[Q.sub.P] = [Q.sub.CM] + [C.sub.QM] x [H.sub.S]. (5)

[Q.sub.CP], [C.sub.QP], [Q.sub.CM], and [C.sub.QM] in (4) and (5) are all values of the previous instant and are known values in the current instant. [Q.sub.CP] and [Q.sub.CM] are determined by the flow and head of the previous instant, respectively, and [C.sub.QP] and [C.sub.QM] are both correlated to the wave velocity and pipe area, respectively, and are positive constants.

The equations for the unit parameters are expressed as

[Q.sub.P] = [Q'.sub.1][D.sup.2.sub.1][square root of [H.sub.P] - [H.sub.S]], (6)

N = [[[n'.sub.1][square root of [H.sub.P] - [H.sub.S]]]/[D.sub.1]], (7)

[M.sub.t] = [M'.sub.1][D.sup.3.sub.1]([H.sub.P] - [H.sub.S)]. (8)

The first derivative of the differential equation of the generator is expressed as

J[d[omega]/dt] = [M.sub.t] - [M.sub.g] - [30[e.sub.g][P.sub.r]/[n.sup.2.sub.r][pi]] [DELTA]n. (9)

Under an off-grid operation, [M.sub.g] = 0 and [e.sub.g] = 0. Thus, the generator equation can be simplified as

n = [n.sub.0] + [[0.1875([M.sub.t] + [M.sub.t0])[DELTA]t]/G[D.sup.2]]. (10)

The guide-vane opening is defined in advance; thus

y = f(t). (11)

Parameter v represents the relative opening; therefore

V = [y/[y.sub.max]]. (12)

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 [H.sub.P] and [H.sub.S] of the inlet and outlet of the rotating wheel, respectively, discharge [Q.sub.P], rotating speed n, torque [M.sub.t], unit speed [n'.sub.1], unit discharge [Q'.sub.1], unit torque [M'.sub.1], guide-vane opening y, and parameters u and v. 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 .

(1) We assume a value for parameter u and obtain the three unit parameters from (3). Inserting X = [square root of [H.sub.P] - [H.sub.S]], [C.sub.1] = [Q.sub.CP]/[C.sub.QP] + [Q.sub.CM]/[C.sub.QM], = 1/[C.sub.QP] + 1/[C.sub.QM], and E = 0.1875[DELTA]t/G[D.sup.2] into (4)-(6) and simplifying the result yield

[X.sup.2] + [C.sub.2] x [Q'.sub.1] x [D.sup.2.sub.1]X - [C.sub.1] = 0. (13)

Then, the two roots are obtained as follows:

[X.sub.1] = [[-[C.sub.2][Q'.sub.1][D.sup.2.sub.1] + [square root of [([C.sub.2][D.sup.2.sub.1][Q'.sub.1]).sup.2] + 4[C.sub.1]]]/2 (14a)

[X.sub.2] = [[-[C.sub.2][Q'.sub.1][D.sup.2.sub.1] - [square root of [([C.sub.2][D.sup.2.sub.1][Q'.sub.1]).sup.2] + 4[C.sub.1]]]/2 (14b)

X represents the square root of the head and is always a positive real number. According to (14a) and (14b), [X.sub.1] > 0 when [C.sub.1] > 0. When [C.sub.1] < 0 and [Q'.sub.1] [less than or equal to] -[square root of -4[C.sub.1]/[C.sup.2.sub.2][D.sup.4.sub.1], both [X.sub.1] and [X.sub.2] are positive real numbers, and we need to select the right one. Under other conditions, when both [X.sub.1] and [X.sub.2] are nonpositive real numbers, Step (1) is repeated.

(2) We insert X into (7) and (8) and obtain rotating speed n and torque [M.sub.t]. We insert [M.sub.t] into (10) and obtain rotating speed [n.sup.*].

(3) We define [OMEGA] = n - [n.sup.*]. When [OMEGA] = 0, parameter u is reasonably assumed. X 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 X into (4)-(6) to obtain the values of [H.sub.P], [H.sub.S], and [Q.sub.P].

The flowchart is shown in Figure 3. When [C.sub.1] > 0, the search interval of [OMEGA] = 0 is the whole equal-opening curve (i.e., from u = 0 to u = 1); when [C.sub.1] < 0, the search interval is the part in which [Q'.sub.1] [less than or equal to] -[square root of -4[C.sub.1]/[C.sup.2.sub.2][D.sup.4.sub.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 [OMEGA] = 0 in a known equal-opening curve and (2) distinguishing between [X.sub.1] and [X.sub.2] when [C.sub.1] < 0 and [Q'.sub.1] [less than or equal to] -[square root of 4[C.sub.1]/[C.sup.2.sub.2][D.sup.4.sub.1]].

Assumption 1. [OMEGA] is a monotonic function in the search segment of the equal-opening curve with parameter u.

Assumption 2. [OMEGA] has opposite signs at the two boundary points of this segment.

From Assumptions 1 and 2, a unique point must exist where [OMEGA] = 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 [X.sub.1] or [X.sub.2], if [OMEGA] is a monotonic function, by inserting [X.sub.1] or [X.sub.2] into [OMEGA], the value that can make the signs of [OMEGA] 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.

3. Search Direction of X

3.1. Directional Derivative of [OMEGA]. According to the directional derivative theorem, the directional derivative of [OMEGA] in line segment [P.sub.j][P.sub.i+1] shown in Figure 2 can be expressed as 

[mathematical expression not reproducible]. (15)

We define [DELTA]L = [square root of [([DELTA][n'.sub.1]).sup.2] + [([DELTA][Q'.sub.1]).sup.2] + [([DELTA][M'.sub.1]).sup.2]] as the length of the segment, cos [alpha] = [DELTA][n'.sub.1]/[DELTA]L, cos [beta] = [DELTA][Q'.sub.1]/[DELTA]L, and cos [gamma] = [DELTA][M'.sub.1]/[DELTA]L as the cosines of the angles separately formed by segment [P.sub.i][P.sub.i+1] and the three coordinate axes in space. Then

[mathematical expression not reproducible]. (16)

Therefore, the sign of [partial derivative][OMEGA]/[partial derivative]u depends on the sign of the sum of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1], [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1], and [DELTA][M'.sub.1][partial derivative][OMEGA]/[partial derivative][M'.sub.1]. When [partial derivative][OMEGA]/[partial derivative]u > 0, [OMEGA] is a monotonically increasing function within [P.sub.i][P.sub.i+1]; when [partial derivative][OMEGA]/[partial derivative]u < 0, [OMEGA] is monotonically decreasing. From the calculation method presented in Section 2.2 and (7)-(10), [OMEGA] can be expressed as

[mathematical expression not reproducible]. (17)

The following parts will consider X = [X.sub.1] as an example for analysis. X = [X.sub.2] can be analyzed by the same method. The partial derivatives of (17) are as follows:

[mathematical expression not reproducible]. (18)

Evidently, [partial derivative][OMEGA]/[partial derivative][n'.sub.1] and [partial derivative][OMEGA]/[partial derivative][M'.sub.1] are constants that are greater than zero, and the sign of [partial derivative][OMEGA]/[partial derivative][Q'.sub.1] varies under the following scenarios. (1) When the operating point is located in the pump region or pump braking regionwhere [n'.sub.1] < 0 and [M'.sub.1] > 0, [partial derivative][OMEGA]/[partial derivative][Q'.sub.1] is positive. (2) When the operating point is located in the turbine braking region or reverse-pump region where [n'.sub.1] > 0 and [M'.sub.1] < 0, [partial derivative][OMEGA]/[partial derivative][Q'.sub.1] is negative. (3) When the operating point is located in the turbine region and [n'.sub.1] < 2E[D.sup.4.sub.1][M'.sub.1]X (the estimated magnitude of [n'.sub.1] is not greater than [10.sup.0]), [partial derivative][OMEGA]/[partial derivative][Q'.sub.1] is positive; otherwise, [partial derivative][OMEGA]/[partial derivative][Q'.sub.1] is negative.

The signs of [DELTA][n'.sub.1], [DELTA][Q'.sub.1], and [DELTA][M'.sub.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  (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 [n'.sub.1] = 2E[D.sup.4.sub.1][M'.sub.1]X 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 [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1], [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1], and [DELTA][M'.sub.1][partial derivative][OMEGA]/[partial derivative][M'.sub.1] along the equal-opening curve.

3.2. Determination of the Monotonicity of [OMEGA]

3.2.1. Comparison between [absolute value of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] and [absolute value of [DELTA][M'.sub.1][partial derivative][OMEGA]/[partial derivative][M'.sub.1]]. Assuming [K.sub.1] as the ratio between [absolute value of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] and [absolute value of [DELTA][M'.sub.1][partial derivative][OMEGA]/[partial derivative][M'.sub.1]], we let

[F.sub.1] = [absolute value of [DELTA][n'.sub.1]/[DELTA][M'.sub.1]] [F.sub.2] = E[D.sup.4.sub.1]X. (19)

According to (18)

[K.sub.1] = [[F.sub.1]/[F.sub.2]]. (20)

The magnitude of E[D.sup.4.sub.1] is usually approximately [10.sup.-4], and the magnitude of X is usually in the range of [10.sup.1] - [10.sup.2]. Hence, the magnitude of E[D.sup.4.sub.1]X is below [10.sup.-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 [F.sub.1] < [F.sub.2], [absolute value of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] < [absolute value of [DELTA][M'.sub.1][partial derivative][OMEGA]/[partial derivative][M'.sub.1]]; however, these areas are also the segments where [absolute value of [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1]] > [absolute value of [DELTA][M'.sub.1][partial derivative][OMEGA]/[partial derivative][M'.sub.1]]. According to the signs of these three terms in Table 1, we know that [partial derivative][OMEGA]/[partial derivative]u > 0. For the other intervals of the equal-opening characteristic curve, [absolute value of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] [much greater than] [absolute value of [DELTA][M'.sub.1][partial derivative][OMEGA]/[partial derivative][M'.sub.1]].

Consequently, except for the BC segment and the neighborhoods of points H and K, the sign of [partial derivative][OMEGA]/[partial derivative]u depends mainly on the signs of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1] and [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1], whereas [DELTA][M'.sub.1][partial derivative][OMEGA]/[partial derivative][M'.sub.1]] can be ignored.

3.2.2. Comparison between [absolute value of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] and [absolute value of [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1]]. Assuming [K.sub.2] as the ratio between [absolute value of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] and [absolute value of [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1]], we let

[mathematical expression not reproducible]. (21)

According to (18)

[K.sub.2] = [[F.sub.3]/[F.sub.4]]. (22)

When [F.sub.3] > [F.sub.4], [absolute value of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] > [absolute value of [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1]].

Figure 6 shows how the signs of [F.sub.3] and [F.sub.4] change along the equal-opening curve. [F.sub.3] 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. [F.sub.4] reaches its maximum value when [Q'.sub.1] = -[square root of -4[C.sub.1]/[C.sup.2.sub.2][D.sup.4.sub.1] and reaches its minimum value near the [n'.sub.1] axis. In the pump braking and turbine operating regions where [Q'.sub.1] smoothly changes (the range of CG or CH, including EG), [absolute value of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] > [absolute value of [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1]]. At the ends of the pump and reverse-pump regions, [absolute value of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] > [absolute value of [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1]] is also possible; however, the two positions are usually beyond the calculation range of the transient process. At the other positions, [absolute value of [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] < [absolute value of [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1]] exists.

When X = [X.sub.1], opposite signs for [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] and [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1]] occur at two positions. The first position is the curve segment EG where [mathematical expression not reproducible]. Then, [partial derivative][OMEGA]/[partial derivative]u > 0. The second position is the curve segment HK where [mathematical expression not reproducible]. Then, [OMEGA][OMEGA]/[partial derivative]u > 0. At other positions, [DELTA][n'.sub.1][partial derivative][OMEGA]/[partial derivative][n'.sub.1]] and [DELTA][Q'.sub.1][partial derivative][OMEGA]/[partial derivative][Q'.sub.1]] are both greater than zero. Therefore, [partial derivative][OMEGA]/[partial derivative]u > 0 along the characteristic curve, and [OMEGA] is a monotonically increasing function.

When X = [X.sub.2], this condition can be demonstrated in the same manner as X = [X.sub.1]. [OMEGA] is a monotonically decreasing function in the pump region as well as in the reverse-pump region (where [Q'.sub.1] [less than or equal to] -[square root of -4[C.sub.1]/[C.sup.2.sub.2][D.sup.4.sub.1]).

In summary, [OMEGA] is a monotonic function in the search segments; thus, Assumption 1 is verified.

3.3. Signs of Q in the Search Boundary

3.3.1. When [C.sub.1] > 0. If the current operating point ([n'.sub.1,u=p], [Q'.sub.1,u=p]) is on the left side of the [n'.sub.1] = 0 axis, [n'.sub.1,u=0] < [n'.sub.1,u=p] < 0, and [Q'.sub.1,u=0] < [Q'.sub.1,u=p] can be obtained. According to (7) and (14a), [[OMEGA].sub.u=0] < 0 if [X.sub.u=0] > [X.sub.u=p] and [n.sub.u=o] < [n.sub.u=p] < 0. For the end point ([n'.sub.1,u=1], [Q'.sub.1,u=1]) in the reverse-pump region, [n.sub.u=1] > 0 and [[OMEGA].sub.u=1] > 0 if [n'.sub.1,u=1] > 0. As a result, when [C.sub.1] > 0, [OMEGA] is a monotonic function within the range of u = [0,1] and [[OMEGA].sub.u=0] x [[OMEGA].sub.u=1] < 0. Therefore, only one point exists where [OMEGA] = 0 along the entire characteristic curve. The same conclusion can also be drawn when the operating point is on the right side of the [n'.sub.1] = 0 axis.

3.3.2. When [C.sub.1] < 0. Similarly, when X = [X.sub.1], [[OMEGA].sub.u=0] < 0 and [[OMEGA].sub.u=1] > 0. When X = [X.sub.2], [[OMEGA].sub.u=0] > 0 and [[OMEGA].sub.u=1] < 0. '

If 0 < a < b < 1, for the boundary point ([n'.sub.1,u=a], [Q'.sub.1,u=a]) where [Q'.sub.1,u=a] = -[square root of -4[C.sub.1]/[C.sup.2.sub.2][D.sup.4.sub.1] in the pump region, we can obtain the result X = [X.sub.1] = [X.sub.2] = [square root of -[C.sub.1]] by inserting [Q'.sub.1,u=1] into (14a) and (14b). Then, inserting X into (17) yields the following result:

[mathematical expression not reproducible]. (23)

At the position where [absolute value of [n'.sub.1]] 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, [[OMEGA].sub.u=a] < 0 when [C.sub.1] < -[n.sup.2.sub.0][D.sup.2.sub.1]/[n'.sup.2.sub.1]; otherwise, [[OMEGA].sub.u=a] > 0. Similarly, at the boundary point ([n'.sub.1u=b], [Q'.sub.1u=b]) where [Q'.sub.1,u=b] = -[square root of -4[C.sub.1]/[C.sup.2.sub.2][D.sup.4.sub.1] in the reverse-pump region, [[OMEGA].sub.u=b] > 0 when [C.sub.1] < -[n.sup.2.sub.0][D.sup.2.sub.1]/[n'.sup.2.sub.1]; otherwise, [[OMEGA].sub.u=a] < 0.

In summary, when [C.sub.1] < -[n.sup.2.sub.0][D.sup.2.sub.1]/[n'.sup.2.sub.1], [[OMEGA].sub.u=0] x [[OMEGA].sub.u=d] < 0 (or [[OMEGA].sub.u=b] x [[OMEGA].sub.u=1] < 0) in the search boundary of the pump region (or reverse-pump region) only if X = [X.sub.2]. When 0 > [C.sub.1] > -[n.sup.2.sub.0][D.sup.2.sub.1]/[n'.sup.2.sub.1], [[OMEGA].sub.u=0] x [[OMEGA].sub.u=a] < 0 (or [[OMEGA].sub.u=b] x [[OMEGA].sub.u=1] < 0) in the search boundary of the pump region (or reverse-pump region) only if X = [X.sub.1].

Based on the systematic discussions of [C.sub.1] > 0 and [C.sub.1] < 0, two conclusions can be drawn. (1) The signs of [OMEGA] at two boundary points in the search segment are opposite; thus, Assumption 2 is verified. (2) When [C.sub.1] < -[n.sup.2.sub.0][D.sup.2.sub.1]/[n'.sup.2.sub.1], X = [X.sub.2]; meanwhile, when [C.sub.1] > -[n.sup.2.sub.0][D.sup.2.sub.1]/[n'.sup.2.sub.1], X = [X.sub.1]. 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. . 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. .

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 Q 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 [OMEGA] = 0, the iteration in the numerical solution was completed.

This study has shown that n 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 n 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 [OMEGA] = 0; 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.

Nomenclature

[A.sub.1], [A.sub.2]: Interpolated coefficients of [n'.sub.1]-[Q'.sub.1] characteristics

[B.sub.1], [B.sub.2]: Interpolated coefficients of [n'.sub.1]-[M'.sub.1] characteristics

[D.sub.1]: Diameter of runner inlet [m]

[e.sub.g]: Self-regulation coefficient of power network

G[D.sup.2]: Flywheel moment [kg x [m.sup.2]]

J: Rotating inertia (= G[D.sup.2]/4g) [kg x m x [m.sup.2]]

[M.sub.g]: Resistance moment of generator [kg x m]

[M.sub.t]: Driving force moment of hydraulic turbine [kg x m]

[M'.sub.1]: Unit torque [kg x m]

[N.sub.QE]: Specific speed (m x kW)

n: Rotational speed [r/min]

[n'.sub.1] : Unit speed [r/min]

[P.sub.r]: Rated power [MW]

Q: Cross-section discharge [[m.sup.3]/s]

[Q'.sub.1]: Unit discharge [[m.sup.3]/s]

t: Time [s]

u: Position of an operating point in opening curve

v: Relative opening degree

X: Square root of head [[m.sup.1/2]]

y: Absolute guide-vane opening [mm] or [[degrees]]

[omega]: Angular velocity (= [pi]n/30) [rad/s].

Subscripts

0: Value of last time step

g: Generator parameter

max: Maximum value

r: Rated value

S: Head of the runner outlet

t : Hydraulic turbine parameter

P: 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.

Acknowledgments

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

References

 J. Yang, K. Zhao, and L. Li, "Analysis of the causes of accident for Unit 7 and 9 in Sayano-Shushenskaya hydropower station," in Proceedings of the Advances in Sciences and Engineering, SREE Conference on Engineering Progress, vol. 2, pp. 17-26, Hong Kong, December 2010.

 C. Nicolet, S. Alligne, A. Bergant, and F. Avellan, "Simulation of water column separation in Francis pump-turbine draft tube," IOP Conference Series: Earth and Environmental Science, vol. 15, no. 2, Article ID 022002, 2012.

 Y. J. Fang and J. Koutnik, "The numerical simulation of the delayed load rejection of a pump-turbine powerplant," IOP Conference Series: Earth and Environmental Science, vol. 15, no. 2, Article ID 022018, 2012.

 W. Zeng, J. Yang, and W. Guo, "Runaway instability of pump-turbines in S-shaped regions considering water compressibility," Journal of Fluids Engineering, vol. 137, no. 5, Article ID 051401, 2015.

 Q. Liu, L. Suo, D. Liu, and H. Chen, Hydropower Station, China WaterPower Press, Beijing, China, 3rd edition, 2007 (Chinese).

 C. Duan, "An analytical formula for calculating water hammer in conduit of hydraulic turbine," A Monthly Journal of Science, vol. 25, no. 3, pp. 248-254, 1980.

 R. W. Angus, "Water hammer in pipes, including those supplied by centrifugal pumps: graphical treatment," Proceedings of the Institution of Mechanical Engineers, vol. 136, no. 1, pp. 245-331, 1937.

 C. Jaeger, Fluid Transients in Hydroelectric Engineering Practice, Blackie, Glasgow, Scotland, 1977.

 E. B. Wylie and V L. Streeter, Fluid Transient, McGraw-Hill International Book, New York, NY, USA, 1978.

 D. E. Goldberg and E. B. Wylie, "Characteristics method using time-line interpolations," Journal of Hydraulic Engineering, vol. 109, no. 5, pp. 670-683, 1983.

 C. Lai, "Comprehensive method of characteristics models for flow simulation," Journal of Hydraulic Engineering, vol. 114, no. 9, pp. 1074-1097, 1988.

 J.-C. Yang and E.-L. Hsu, "Time-line interpolation for solution of the dispersion equation," Journal of Hydraulic Research, vol. 28, no. 4, pp. 503-520, 1990.

 B. W. Karney and M. S. Ghidaoui, "Flexible discretization algorithm for fixed-grid MOC in pipelines," Journal of Hydraulic Engineering, vol. 123, no. 11, pp. 1004-1011, 1997.

 M. H. Chaudhry, Applied Hydraulic Transients, Springer, New York, NY, USA, 3th edition, 2014.

 M. H. Chaudhry and M. Y. Hussaini, "Second-order accurate explicit finite-difference schemes for waterhammer analysis," Journal of Fluids Engineering, vol. 107, no. 4, pp. 523-529, 1985.

 V. Guinot, "Riemann solvers for water hammer simulations by Godunov method," International Journal for Numerical Methods in Engineering, vol. 49, no. 7, pp. 851-870, 2000.

 Y.-H. Hwang and N.-M. Chung, "A fast Godunov method for the water-hammer problem," International Journal for Numerical Methods in Fluids, vol. 40, no. 6, pp. 799-819, 2002.

 M. Zhao and M. S. Ghidaoui, "Godunov-type solutions for water hammer flows," Journal of Hydraulic Engineering, vol. 130, no. 4, pp. 341-348, 2004.

 C. Wang and J.-D. Yang, "Water hammer simulation using explicit-implicit coupling methods," Journal of Hydraulic Engineering, vol. 141, no. 4, Article ID 04014086, 2014.

 J. Zhan and Z. He, "Mathematical model and simulation of large hydro turbine," Water Power, vol. 31, no. 12, pp. 54-55, 2005 (Chinese).

 J. Chang, Transition Process of Hydraulic Machinery, Higher Education Press, Beijing, China, 2005 (Chinese).

 S. Derakhshan and A. Nourbakhsh, "Theoretical, numerical and experimental investigation of centrifugal pumps in reverse operation," Experimental Thermal and Fluid Science, vol. 32, no. 8, pp. 1620-1627, 2008.

 X.-Q. Li, J.-S. Chang, C.-S. Li, P. Chen, and X.-L. Tang, "Hydraulic transient controlling rules coupled with ball-valve and guide-vane in pumped storage power station," Journal of Hydraulic Engineering, vol. 45, no. 2, pp. 220-226, 2014 (Chinese).

 A. P. Boldy, "Representation of characteristics of reversible pump turbines for use in waterhammer simulation," in Proceedings of the 4th International Conference on Pressure Surges, pp. 287-296, Bath, UK, September 1983.

 N. Walmsley, The numerical representation of pump-turbine performance characteristics [Ph.D. thesis], University of Warwick, Warwicks, UK, 1986.

 W. Weng and J. Yang, "Study on the partition of the characteristic curve and the numerical simulation of the boundary of the unit," Water Resources and Hydropower Engineering, vol. 34, no. 2, pp. 54-56, 2003 (Chinese).

 J. Yang, J. Yang, and C. Wang, "Mathematical model and simulation of pump turbine with characteristic space curves," Journal of Hydroelectric Engineering, vol. 32, no. 5, pp. 244-250, 2013 (Chinese).

 J. Yang and J. Yang, "B-spline surface construction for the complete characteristics of pump-turbine," in Proceedings of the Asia-Pacific Power and Energy Engineering Conference (APPEEC '12), pp. 1-5, Shanghai, China, March 2012.

 D. Varberg, E. J. Purcell, and S. E. Rigdon, Calculus, Machinery Industry Press, Beijing, China, 1st edition, 2009.

 J. Liu, Study on the variety of characteristic curve and its effect on hydraulic transient process in pumped storage plant [M.S. thesis], Wuhan University, Wuhan, China, 2005 (Chinese).

 W. Zeng, J. Yang, J. Hu, and J. Yang, "Guide-vane closing schemes for pump-turbines based on transient characteristics in S-shaped region," Journal of Fluids Engineering, vol. 138, no. 5, Article ID 051302, 2016.

Yanna Liu, Jiandong Yang, Jiebin Yang, Chao Wang, and Wei Zeng

State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China

Correspondence should be addressed to Jiebin Yang; jbyang830212@whu.edu.cn

Received 11 March 2016; Accepted 31 May 2016

http://dx.doi.org/10.1155/2016/1504659

Caption: Figure 1: Characteristic grid.

Caption: Figure 2: Nonuniform B-spline surface of the turbine characteristics.

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

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

Caption: Figure 5: Variation in [F.sub.1]([n'.sub.1], [Q'.sub.1]).

Caption: Figure 6: Variation in [F.sub.3]([n'.sub.1], [Q'.sub.1]) and [F.sub.4]([n'.sub.1], [Q'.sub.1]).

Caption: Figure 7: Schematic diagram of the piping system of a hydropower station.

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

Caption: Figure 9: Pumped-storage station model.

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

Caption: Figure 11: Comparison of the simulation and field test.
```Table 1: Signs of the different terms at different segments.

[[partial derivative][OMEGA]/
Segment      [[partial derivative]      [DELTA][n'.sub.1]
[n'.sub.1]]

Pump region     AB                   +                        +
BC
+                        +
Pump braking    CD                   +                        +
region

Turbine         DE                   +                        +
region        EF                   +                        +
FG                   +                        +
GH                   +                        +
HI                   +                        -

Braking         IJ                   +                        -
region

Reverse-pump    JK                   +
region        KL                   +                        +

[[partial derivative][OMEGA]/
Segment  [[partial derivative][n'.sub.1]]
[DELTA][n'.sub.1]

Pump region     AB                     +
BC
+
Pump braking    CD                     +
region

Turbine         DE                     +
region        EF                     +
FG                     +
GH                     +
HI                     -

Braking         IJ                     -
region

Reverse-pump    JK
region        KL                     +

[[partial derivative][OMEGA]/
Segment      [[partial derivative]      [DELTA][Q'.sub.1]
[Q'.sub.1]]

Pump region     AB                   +                        +
BC
+                        +
Pump braking    CD                   +                        +
region

Turbine         DE                   +                        +
region        EF                   -                        +
FG                   -                        +
GH                   -                        -
HI                   -                        -

Braking         IJ                   -                        -
region

Reverse-pump    JK
region        KL                   -                        -

[[partial derivative][OMEGA]/
Segment  [[partial derivative][Q'.sub.1]]
[DELTA][Q'.sub.1]

Pump region     AB                     +
BC
+
Pump braking    CD                     +
region

Turbine         DE                     +
region        EF                     -
FG                     -
GH                     +
HI                     +

Braking         IJ                     +
region

Reverse-pump    JK                     +
region        KL                     +

[[partial derivative][OMEGA]/
Segment      [[partial derivative]      [DELTA][M'.sub.1]
[M'.sub.1]]

Pump region     AB                   +                        -
BC
+                        -
Pump braking    CD                   +                        +
region

Turbine         DE                   +                        +
region        EF                   +                        +
FG                   +                        -
GH                   +                        -
HI                   +                        -

Braking         IJ                   +                        -
region

Reverse-pump    JK                   +
region        KL                   +                        -

[[partial derivative][OMEGA]/
Segment  [[partial derivative][M'.sub.1]]
[DELTA][M'.sub.1]

Pump region     AB                     -
BC
-
Pump braking    CD                     +
region

Turbine         DE                     +
region        EF                     +
FG                     -
GH                     -
HI                     -

Braking         IJ                     -
region

Reverse-pump    JK
region        KL                     -

Table 2: Basic pump-turbine information of the engineering cases.

Rated P   Rated n   Rated H      Rated Q      Rotational inertia
(MW)       (rpm)      (m)     ([m.sup.3]/s)    (kg x [m.sup.2])

306.1      428.6      419         83.64            6000000

Rated P   Calculated     Calculated      Calculated
(MW)        H (m)      Q ([m.sup.3]/s)     P (MW)

306.1       457.00          73.18          306.1

Table 3: Basic unit information of the pipe system.

Pipe section              Inlet   Upstream     Upstream
main pipe   branch pipe

Equivalent diameter (m)   0.597     0.357        0.202
Length (m)                2.251    24.154        2.195

Pipe section               Upstream     Downstream   Outlet
branch pipe   main pipe

Equivalent diameter (m)      0.300        0.426      0.675
Length (m)                   4.891        16.285     1.062

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

(m-kW)     Inlet      Outlet    Guide-vane   Number of    Number of
diameter   diameter     height      blades     guide vanes
(mm)       (mm)        (mm)

37.91       280       146.34      24.44          9           20

(m-kW)    Rated n   Rated H   Rated Q     Rotational
(rpm)      (m)      (L/s)        inertia
(N x [m.sup.2])

37.91      1000      10.54     49.1          66.4

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

Rated P   Rated n   Rated H      Rated Q      Rotational inertia
(MW)       (rpm)      (m)     ([m.sup.3]/s)    (kg x [m.sup.2])

306.1       500      517.4        66.72            6000000

Rated P   Tested H     Tested Q      Tested P (MW)
(MW)        (m)      ([m.sup.3]/s)

306.1      535.03        62.78            298
```
Title Annotation: Printer friendly Cite/link Email Feedback Research Article Liu, Yanna; Yang, Jiandong; Yang, Jiebin; Wang, Chao; Zeng, Wei Mathematical Problems in Engineering Report 9CHIN Jan 1, 2016 7713 A Chaos Robustness Criterion for 2D Piecewise Smooth Map with Applications in Pseudorandom Number Generator and Image Encryption with Avalanche... Service and Price Decisions of a Supply Chain with Optional After-Sale Service. Boundary layer Boundary layers (Fluid dynamics) Hydraulic turbines Hydroelectric power Hydroelectric power plants Turbines Water-power