Stability analysis of governor-turbine-hydraulic system by state space method and graph theory

Energy 114 (2016) 613e622

Contents lists available at ScienceDirect

Energy journal homepage: www.elsevier.com/locate/energy

Stability analysis of governor-turbine-hydraulic system by state space method and graph theory Xiaodong Yu*, Jian Zhang, Chengyu Fan, Sheng Chen College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing, People's Republic of China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 13 March 2015 Received in revised form 29 July 2016 Accepted 30 July 2016

The coefficient matrix of the state equations is essential for stability analysis of the governor-turbinehydraulic (GTH) system by the state space method. With plant layout becoming more and more complicated, it is important to derive the state matrix quickly and accurately. Based on the stability analysis theory of the GTH system, this paper investigates regular features of the state equations that describe small fluctuations in the state variables of the system. The equations for unsteady flow in the pipeline system are conveniently given by using graph theory. By specifying the order of the state variables and using matrix transformation, an innovative method for solving the coefficient matrix of the state equations is established, and the stable region of the system can be given with the eigenvalue method. The proposed method is used to analyze the stability of a practical hydropower station during small fluctuation, which is also verified in the hydraulic transient model of hydropower system on the basis of characteristics method. © 2016 Elsevier Ltd. All rights reserved.

Keywords: Hydropower station Stability analysis State matrix Small fluctuation Modeling

1. Introduction Hydropower is a major renewable energy resource, which is widely utilized all over the world. With the rapid exploitation of hydropower resources in China, a lot of hydropower stations are constructed and the gross installed hydropower capacity of China ranked first in the world [1]. Stable operation of hydro-units of a large scale hydropower station has great influence on the safety of units [2] and even the whole power system [3]. So it is an important issue on the design of a hydropower generation system to ensure its stability and regulating quality. It is generally considered that the GTH system of a hydropower station is steady when the units are connected to a large power grid under a small disturbance condition. However, small and mini hydropower stations are mostly isolated from national grid, which make a significant contribution to energy needs [4]. In addition, with the application of the highvoltage direct current (HVDC) transmission technology, large power stations may be in an islanding condition as well. Particularly, some hydropower stations in west China are diversion-type stations with a fairly long headrace tunnel and a surge tank, and the

* Corresponding author. E-mail addresses: [email protected] (X. Yu), [email protected] (J. Zhang), [email protected] (C. Fan), [email protected] (S. Chen). http://dx.doi.org/10.1016/j.energy.2016.07.164 0360-5442/© 2016 Elsevier Ltd. All rights reserved.

stable issue is more prominent [5]. Regarding these hydropower stations, the GTH system should be designed to be stable in isolated operation during design process. Governor-turbine-hydraulic system is complex because of strong couplings of hydraulic, mechanical and electrical system, which is an important application area for control engineering [6]. The schematic of hydropower plant and GTH system are illustrated in Fig. 1. Fang [7] introduced the basic models for hydraulic turbine and control system. The speed oscillations following a load change are stable or unstable depending upon the values of the parameters of waterway system, hydro-unit and governor, so that, the sensitivity analysis of such parameters would help to optimize system design or the turning of governors. There are in general two approaches to analyze the stability of GTH system. One is numerical simulation in time domain [8]. Although the numerical simulation method can be taken as many nonlinear factors of the system as possible, the mathematical model is complex and programming is also very difficult. The results are affected by numerical errors, particularly in a small disturbance condition, which might result in wrong judgments. In addition, it is not convenient to optimize system parameters or the turning of governors. The second one uses the automatic control theory. The state space method and the transfer function method have been widely used for studying the stability and performance of the GTH system.

614

X. Yu et al. / Energy 114 (2016) 613e622

upstream surge tank

headrace tunnel Unit

downstream

penstock tailrace tunnel

(a) General layout Plant Hydraulic System

Turbine Generator

Grid

Speed deviation

Wicket Gate

Governor

Speed setting

(b) GTH system Fig. 1. Schematic of hydropower plant (a) general layout; (b) GTH system.

Chaudhry [9] used the Routh-Hurwitz stability criterion to analyze the stability of the GTH system. Hagihara et al. [10] established the stable region of turbine-generator units with proportional integral differential (PID) governor by using the root locus method. Thorne et al. [11] applied the state space method and eigenvalue analysis to examine the effects of system parameters on the stability boundaries of single hydro unit connected to a large power system. Kishor [12] used the state space hydro plant model with inelastic and elastic water column effect to simulate the regulation of turbine speed with the linear quadratic approach. Chen et al. [13] introduced a nonlinear model based on state equations to study the nonlinear behaviors of GTH system and a fuzzy sliding mode governor was designed. Load rejection transients were simulated with a nonlinear model of a hydro-turbine governing system in [14], and the stable regions were presented by analyzing the coefficient matrix of the system. Zhou et al. [15] utilized a higher-order elastic model of pipe flow in state equations. Xu et al. [16] introduced fractional order calculus to the state space model of GTH system and the stability was studied. Guo et al. [17] studied the regulation quality for frequency response of turbine regulating system of an isolated hydroelectric power plant, and the expressions of coefficients in overall transfer function was derived. Yang et al. [18] proposed a linear state space model of the hydro-turbine governing system with an open tailrace channel, and Guo et al. [19] used the Hopf bifurcation theory to study the stability of the system. Donaisky et al. [20] studied the hydraulic amplifier of a governor and the frequency control performance of the hydropower plant is improved. Martínez-Lucas et al. [21] investigated the performance of the frequency control of isolated power systems including a hydropower plant. Based on the transfer function model of the GTH system, a lot of methods for improving the setting performance of the system parameters were presented, such as the robust control methods [22], improved

particle swarm optimization algorithm [23], adaptive grid particle swam optimization [24], improved gravitational search algorithm [25], fuzzy sliding mode control [26]. However, most of the GTH system used in the literature are simple, including a turbine-generator unit, a single penstock and an upstream reservoir. Even though the system is simple, the authors spent a lot of time to derive the coefficient expressions of the state space model, which is very complex and easy to make mistake. But in a practical hydropower plant, the GTH system is much more complex (i.e. multi-unit, hydraulic coupling and surge tanks), let along multi-hydropower plants connection. When the system is complicated, both the state space model and the transfer function model are not easy to establish and it is also difficult to judge the stability of a nonlinear and multi-variable system with the automatic control theory directly. The state space method is used to analyze the stability of the GTH system under small disturbance in this paper. The coefficient matrix of the state equations that describe the small fluctuations in the state (dependent) variables of the system is essential for stability analysis of the GTH system. As we have discussed before, how to get the coefficient matrix of the state equations accurately and quickly is significant in stability analysis of a rather complicated GTH system with fewer simplifications. Inspired by the above discussions, this paper investigates the regular features of the state equations that describe the small fluctuations in the GTH system, and the equations for unsteady flow in the complex pipeline system are conveniently given by using graph theory. By specifying the order of the state variables and using matrix transformation method, an innovative method for solving the coefficient matrix of the state equations is established. With this method, the stability region of the system can be given by using the eigenvalue method, which can be used for optimizing the system parameters or the turning of governors. Then, the proposed

X. Yu et al. / Energy 114 (2016) 613e622

method is used to analyze the stability of a practical hydropower station during small fluctuations, which is also verified in the hydraulic transient model of hydropower system on the basis of characteristics method. 2. Governing equations The following assumptions are made to develop the governing equations: 1. The steady-state characteristics of turbine hold for unsteadystate situations and homologous relationships are valid. 2. The changes in the turbine speed, head and gate opening are small; thus nonlinear relationships can be assumed linear. 3. Rigid water column theory and quasi-steady friction terms are applied. 4. The systems are under islanding conditions. 5. The governor has no dead band, backlash or hysteresis. By making these assumptions, the following equations may be written for the components of the GTH system in a hydropower station, as shown in Fig. 1. 2.1. Turbine and generator equations The flow passing through the turbine and the power produced by the turbine are given by

pffiffiffiffi Q ¼ D21 Q11 H P ¼ T11 D1 H 

(1) np 30

(2)

where Q ¼ the discharge, H ¼ the head across the turbine, P ¼ the turbine output, n ¼ the rotating speed, D1 ¼ the runner diameter, and T ¼ the torque. Subscript 11 refers to the unit value. Both Q11 and T11 are functions of the unit speed, n11, and the wicket gate position, t, as given by the characteristic curves of the turbine. By using Taylor series expansion, Eqs. (1) and (2) can be linearized in dimensionless form:

Tm

615

d4 ¼ p  pG ¼ p  x  sp 4 dt ½GD2 n2

in which Tm ¼ 365P00 ¼ the mechanical starting time with [GD2] ¼ the moment of inertia of rotating fluid and mechanical parts in the turbine-generator unit in tons-m2, n ¼ the rotating speed in rpm, and P ¼ the power in kilowatts. The power absorbed by the generator, PG, depends on the change in the external load, X, vpG G0 0 , x ¼ XX and on the speed variation. Defining pG ¼ PGPP X0 ,sp ¼ v4 G0

gives pG¼xþsp4. Taking Eqs. (5) and (7) and eliminating p, the turbine-generator equation is expressed as

    s9i  spi d4i s s 1 ¼ x 4i þ 10i mi þ 8i qi þ  Tmi i dt Tmi Tmi Tmi

2.2. Governor equation The classical PID is the most common form of controller used in governing. The PID governor equation, which relates the speed change to the position of the main control mechanism, may be expressed as

  dm d2 4 d4 4 bt þ bp Td þ bp m ¼ Td Tn 2  ðTn þ Td Þ dt dt dt

"

 #   s9i  spi d4i þ bti þ bpi Tdi þ Tni Tni þ Tdi þ Tni  Tdi  dt Tmi

s dmi s dqi þ Tni  Tdi  8i  Tdi  10i Tmi dt Tmi dt ¼ 4i  bpi mi (10)

in which i ¼ 1,2, …,N, N ¼ number of the turbine-generators.

h s i p ¼ ð1 þ s3 Þ4 þ s4 m þ 1  3 x 2

(4)

2.3. Unsteady flow equation in pipeline system

* vT11 vn*11

* vT11 v *11

* vQ11 , vn*11 * t ¼

* ¼ Q11 , T * ¼ T11 , n* ¼ n11 , s2 ¼ t , s3 ¼ , s4 ¼ t with Q11 n110 11 T110 Q110 11 t are parameters that depend on the turbine characteristic curves t0

and the initial operating point, which can be derived by interpolation and finite difference (FD) method. Subscript 0 refers to initial 2 , 1 2 steady state. Defining s5 ¼ 1s s6 ¼ 2s s7 ¼ 2s 1s1 , 1s1 ,  1   s9 ¼ 1  s23  s6 þ 1 þ s3 , and s8 ¼ 1  s23  s5 ,   s10 ¼ 1  s23  s7 þ s4 , Eqs. (3) and (4) can be written as

p ¼ s8 q þ s9 4 þ s10 m

(5)

x ¼ s5 q þ s6 4 þ s7 m

(6)

The dimensionless form of the torque equation is

(9)

in which bp ¼ the permanent speed drop of the governor; bt ¼ its temporary speed drop; Td ¼ its dashpot time constant; and Tn ¼ its promptitude time constant. 2 Differentiating Eq. (8) with respect to t, and eliminating ddt42 , the PID governor equation is expressed as

(3)

* vQ11 v *11

(8)

in which i ¼ 1,2, …,N, N ¼ number of the turbine-generators.

1 q ¼ s1 4 þ s2 m þ ð1  s1 Þx 2

Q Q0 tt0 PP0 HH0 0 in which4 ¼ nn n0 ,m ¼ t0 ,q ¼ Q0 ,p ¼ P0 ,x ¼ H0 , and s1 ¼

(7)

As the changes in the dependent variables are limited to small values during small fluctuations, rigid water column theory and quasi-steady friction terms are applied to analyze the unsteady flow in pipeline. The unsteady flow for a simple pipeline is described by

L dQ ¼ Hu  aQ 2  Hd gA dt

(11)

in which L ¼ the length of the pipe, A ¼ the cross-sectional area of the pipe, a ¼ the head loss coefficient alone the pipe, Q ¼ the flow in the pipe, Hu ¼ the upstream piezometric head; and Hd ¼ the downstream piezometric head. Different hydropower stations have different layouts, so it is difficult to give a general formula for pipeline system as the turbine-generator Eq. (8) or the governor Eq. (10). In this section, the regular features of the unsteady flow equation are investigated,

616

X. Yu et al. / Energy 114 (2016) 613e622

and the equations for unsteady flow in the pipeline system are conveniently given by using graph theory. As shown in Fig. 2, every four units have been installed in one water conservancy system with 11 pipe sections. According to the principle of continuity, the state variables regarding discharge are marked in each pipe in the diagram. It is obviously that there are 12 state variables for analyzing the stability of the GTH system, so 12 first order differential equations are needed to deduce the state equations of the system. There are four units in the system, so we can derive 8 equations from the turbine-generator Eq. (8) and the governor Eq. (10), while the other four differential equations regarding state variables are still needed. Take 1# unit's water diversion system as example. By coupling the dynamic Eq. (11) of the pipes (i.e. 1#, 2#, 3# and 4#) between the upstream reservoir and the downstream reservoir, and simplifying and linearizing these equations at the initial steady state, we can derive the following equation (Due to space limitations, the specific steps are not presented.)

ðB1 þ B2 þ B3 þ B4 ÞC1;1

state variables, all the 12 state equations the GTH system shown in Fig. 2 can be solved out.

2.4. Differential equations of other state variables Besides the dimensionless rotating speed 4, the dimensionless wicket gate opening m and dimensionless flow rate q, there are other state variables in stability analysis, such as the water level of the surge tank and its equation is

dZsu Asu ¼ Qsu dt

(14)

in which, Zsu ¼ the water level of the surge tank, Asu ¼ the crosssectional area of the surge tank and Qsu ¼ the water flow into the surge tank.

dq1 dq dq dq þ ðB1 þ B2 ÞC1;2 2 þ B1 C1;3 3 þ B1 C1;4 4 ¼ðD1 þ D2 þ D3 þ D4 ÞC1;1 q1 þ ðD1 þ D2 ÞC1;2 q2 þ D1 C1;3 q3 dt dt dt dt þ D1 C1;4 q4  x1 (12)

Li in which, Bi ¼ gA , i ¼ 1~K(K is the total number of the pipeline i

3. Stability criteria and coefficient matrix of state equations

Q

sections), Ci;j ¼ Hj0i0 , i ¼ 1~N, j ¼ 1~NQ9(N is number of the units and 0

NQ is number of the flow variables) and Di ¼ 2ai Qi0 , i ¼ 1~K, Subscript 0 refers to initial steady state value. Comparing Eq. (12) and Fig. 2, it can be found that the coefficients of the state variables are corresponding to the layout of pipes. As shown in Fig. 2, the state variable q1 is corresponding to the pipe Section 1#, 2#, 3# and 4#; q2 is corresponding to the pipe Section 1# and 2#; q3 is corresponding to the pipe Section 1#, and q4 is corresponding to the pipe Section 1#. Therefore, the dynamic equation of pipes corresponding to 2# unit can be derived directly according to Fig. 2.

ðB1 þ B2 ÞC2;1

The stability of the GTH system of a hydropower station depends on the eigenvalues of coefficient matrix A (ln¼snþiun, where s ¼ the attenuation factor, u ¼ the natural frequencies of free vibrations in the system and n ¼ number of the state variables). If all the real part of eigenvalues are negative, i.e.sn<0, the system is stable, otherwise unstable. The GTH system is in a critical stable state if smax¼0. The eigenvalues are found easily by means of a computer with standard subroutines using QR transformations, so the coefficient matrix A of the state equations is essential for stability analysis.

dq1 dq dq dq þ ðB1 þ B2 þ B5 þ B6 ÞC2;2 2 þ B1 C2;3 3 þ B1 C2;4 4 ¼ðD1 þ D2 ÞC2;1 q1 þ ðD1 þ D2 þ D5 þ D6 ÞC2;2 q2 dt dt dt dt þ D1 C2;3 q3 þ D1 C2;4 q4  x2

Similarly, the dynamic equations of 3# and 4# unit can also be deduced conveniently. By substituting Eq. (6) into dynamic equations and replacing the dimensionless head variation x of unit with q1 q2 q3 q4

1

q3 q4 7 q1 2

q2

q4 10 q3 8 q2 5 q1 3

Fig. 2. Pipeline system of GTH system.

4# 3# 2# 1#

q4 11 q3 9 q2 6 q1 4

(13)

Combined with the governing equations above, the state equations for the GTH system are yielded and may be expressed in matrix form:

dY ¼ AY þ Bx dt

(15)

in which Y ¼ the state vector, A ¼ the coefficient matrix of state equations, B ¼ input matrix, and x is the relative disturbance. Unless the numbers of the state variables of the system are few, the coefficient matrix A can be deduced directly, and it is difficult to solve out the state equations. By left multiplying a coefficient matrix E on both sides of Eq. (15), Eq. (15) can be transformed as:

X. Yu et al. / Energy 114 (2016) 613e622

617

Fig. 3. Flowchart of the present method.

dY ¼ FY þ R E dt

unit's rotational speed variation, the dimensionless unit's wicket gate opening variation, the dimensionless discharge variation and the dimensionless variation of other variables (e.g. the water level in surge tank) respectively. By introducing N ¼ number of the units, NQ ¼ number of the discharge variables (including the unit's discharge), and W ¼ number of the other variables, the total number of the state variables is 2 N þ NQ þ W. Defining M ¼ 2 N þ NQ þ W, M differential equations regarding state variables are needed to analyze the stability of the GTH system. Taking Eqs. ((8), (10), (12) and (14), the governing equations can be written in matrix form as Eq. (16):

(16)

in which, E, F and R are coefficient matrixes of the system (F¼EA and R¼EBx). Compared with matrix A, it is easier to give matrix E, F and R. These three matrixes are regular and easier for programming. When matrix E, F and R are solved out, equation (15) can be deduced directly by left multiplying by the inverse matrix of the coefficient E on both sides of Eq. (16). The state variables 4,m,q and z are specified in order, so as to find out the regular features of matrix E, F and R and program more easily. Those variables 4,m,q and z represent the dimensionless

0

e1;1 … e1;M B « 1 « B B e2N;1 / e2N;M B …………………… Be B 2Nþ1;1 … e2Nþ1;M B « 1 « B B e / e3N;M B …………………… 3N;1 B B e3Nþ1;1 … e3Nþ1;M B @ « 1 « / eM;M eM;1

surge tank

1 q1

3

q2

1# 2#

q1 2 q2 4

q3 5

Fig. 4. Schematic diagram of the hydropower station in case study.

0

Table 1 Pipe system data. Pipe

L (m)

D (m)

a (m/s)

N

n

1# 2# 3# 4# 5#

338.7 132.1 359.5 132.1 395.1

9.0 10.7 9.0 10.7 16.6

1000 1000 1000 1000 900

96 40 102 40 119

0.014 0.012 0.014 0.012 0.014

0 dy1 1 dt C 1B B C B C C CB B C CB C CB C CB C CB C CB C CB CB « C C CB C CB C CB C CB C CB C AB C B C B C @ A dym dt 10 y1 1

f1;1 … f1;M B « 1 « B B f2N;1 / f2N;M B …………………… Bf B 2Nþ1;1 … f2Nþ1;M B « 1 « ¼B B / f3N;M f3N;1 B …………………… B B f3Nþ1;1 … f3Nþ1;M B @ « 1 « / fM;M fM;1

Note: n is Manning roughness coefficient, D is equivalent diameter and N is number of reaches for each pipe.

0 r1 1

CB C B C CB C B C CB C B C CB C B C CB C B C CB C B C CB C B C CB « C þ B « C CB C B C CB C B C CB C B C CB C B C CB C B C A@ A @ A yM

(17)

rM

Table 2 Parameters of surge tank and unit. Surge tank 2

Asu (m ) 1200

Unit 2

Aorf (m ) 62

zin

zout

0.6

0.8

PR (MW) 252.6

HR (m) 81

NR (r/min) 115.4

QR (m3/s) 344.7

WD2/g (t$m2) 60000

D (m) 5.9

Table 3 Parameters of initial operating condition. Turbine

H (m)

Q (m3/s)

t ( )

n11 (rpm)

Q11 (m3/s)

T11 (N m)

s1

s2

s3

s4

1# 2#

81.0 81.0

344.6 344.6

33.5 33.5

75.6 75.6

1.1 1.1

1256.4 1256.1

0.07 0.07

0.74 0.74

0.81 0.81

0.55 0.55

618

X. Yu et al. / Energy 114 (2016) 613e622

Fig. 5. Characteristic curves of the turbine (a) unit discharge; (b) unit torque.

The first part coefficient expressions ei,i¼1, fi;i ¼ fi;iþN ¼

ðs10i Þ Tmi ,

fi;iþ2N ¼

ðs8i Þ Tmi ,

ri ¼

xi Tmi ,

ðs9i spi Þ Tmi ,

eiþN,i¼TniþTdiþTniþTdiþfi,i,

eiþN,iþN¼(btiþbpi)TdiþTniTdifi,iþN, eiþN,iþ2N¼TniTdifi,iþ2N, fiþN,i¼1and fiþN,iþN¼bpi(i¼1~N) are corresponding to the turbine-generator equation (8) and the governor equation (10). When the governor type is PID or PI (Tn ¼ 0) regulation mode, the coefficient expressions are the same. That is, the coefficient expressions in coefficient matrix of this part are the same. All the coefficients that are not assigned are 0, and similarly hereinafter. The second part coefficient expressions e2Nþi,j¼BXCY, f2Nþi,i¼S6,i,

0

f2Nþi,Nþi¼S7,i, f2Nþi,j¼DXCY, and f2Nþi;2Nþi ¼ f2Nþi;2Nþi  S5;i (i ¼ 1~N, j ¼ 2 N þ 1~NQ) are corresponding to the unsteady flow equation in pipeline system on the basis of graph theory, in which the expressions of B,C and D have been introduced before. It should be noted that as the dimensionless unit's head variation x is expressed by the state variables in dynamic equation, the coefficients 0 f2Nþi,i¼S6,i,f2Nþi,Nþi¼S7,i,f2Nþi;2Nþi ¼ f2Nþi;2Nþi  S5;i are intro0 duced, where f2Nþi;2Nþi have been known with the assignment of f2Nþi,j¼DXCY. Now, all the 12 state equations corresponding to Fig. 2 are given with standard coefficient expressions. The remaining part in the coefficient matrix are the coefficients

X. Yu et al. / Energy 114 (2016) 613e622

619

16

12

B (0.5,12)

Td

Unstable

Stable 8

A (0.3,8)

4 0. 2

0. 3

0. 4

0. 5

0. 6

bt Fig. 6. Stable region in bt-Td plane for fixed Tn ¼ 1.5.

Table 4 Eigenvalues with governor parameters A and B. Parameters (bt, Td)

Eigenvalues ln ¼ sn þ iun

A (0.3,8)

s u s u

B (0.5,12)

0.2083 0.4310 0.5110 0.0000

0.2083 0.4310 0.5460 0.0000

0.1734 0.4602 0.0020 0.0669

0.1734 0.4602 0.0020 0.0669

0.3614 0.0000 0.0605 0.1389

0.3179 0.0000 0.0605 0.1389

0.0004 0.0676 0.0747 0.1310

0.0004 0.0676 0.0747 0.1310

Dimentionless speed

103% Numerical simulation State equations

102%

101%

100%

99%

98% 0

100

200

300

400

500

Time (s)

Water level in surge tank (m

1821.4

Numerical simulation State equations

1821.2

1821.0

1820.8

1820.6

1820.4 0

100

200

300

400

500

Time (s) Fig. 7. Variation of dimensionless rotating speed and water level in surge tank with two different models, using governor parameters A¼(0.3,8).

corresponding to state variables despite (4,m,q) of the unit. If there are other discharge variables, the unsteady flow equation can also be directly written out by using graph theory, while the equations

for other state variables can be derived according to the physical meanings (e.g. Eq. (14)). The order of the state variables should be the same as above.

620

X. Yu et al. / Energy 114 (2016) 613e622 103% Numerical simulation State equations

Dimentionless speed

102%

101%

100%

99%

98% 0

100

200

300

400

500

Time (s)

Water level in surge tank (m

1822.4

Numerical simulation State equations

1821.8

1821.2

1820.6

1820.0

1819.4 0

100

200

300

400

500

Time (s) Fig. 8. Variation of dimensionless rotating speed and water level in surge tank with two different models, using governor parameters B¼(0.5,12).

All the coefficients in matrixes E, F and R in Eq. (17) have already been assigned in standard form. When Eq. (17) is left multiplied by the inverse matrix of the coefficient E on both sides, the essential coefficient matrix A of the state equations is solved out. The procedure of the present method is schematically shown in Fig. 3. It should be noted that this part focuses on introducing the method to establish the state space model of the complex GTH system and analyze the stability accurately and quickly, the effect of the synchronous generator and the power grid is described simplistically by Eq. (7). The external load variation is a boundary condition for the GTH system. If we need to consider these factors, such as excitation system and electrical system which are exist in practical hydropower plants, we can modify Eq. (7) and add some differential equations describing the new state variables in Section 2.4. The corresponding coefficient expressions in Eq. (17) are replaced and added, while other parts keep the same, so the propose method is also applicable. 4. Case study The model is used to analyze the stability of the GTH system in a practical hydropower station in China. As shown in Fig. 4, every two-units share a common waterway system, which consists of the intake, the gate shaft, the headrace tunnel, penstock, the tailrace tunnel, the downstream surge tank and so on. The system parameters are presented in Tables 1 and 2. According to Fig. 4, the state variables of the GTH system of this hydropower station are defined as 41,42,m1,m2,q1,q2,qt1 and zsu. By

using the method above, the coefficient matrix A of the state equations is derived conveniently, which is used to judge the stability of the GTH system. Due to space limitations, the derivation process is not presented. In order to investigate the stability of the GTH system, the critical operating condition is used where the upstream water level is 1903.3 m, the downstream water level is 1820.0 m, the units operate at rated condition and the external load disturbance is 2% rated value. The initial parameters of the turbine are known as shown in Table 3 and Fig. 5. Substituting these parameters into Eq. (17) and left multiplying the inverse matrix of the coefficient E on both sides, the coefficient matrix A of the state equations is solved out. The stable region of the system can be given with the eigenvalue method, which reveals the influence of the parameters on system stability. The stable region in bt-Td plane for fixed Tn ¼ 1.5 are shown in Fig. 6. To verify the proposed method, two models are used to simulate the small fluctuation condition in the GTH system. In the first model, the state equations are solved out by using the fourth-order RungeeKutta method. So the time history of the absolute value of the state variable during the small disturbance is given. The second one is the hydraulic transient model of hydropower system on the basis of characteristics method. This model has been introduced in Ref. [27] in detail, which is verified by the field test as well. The computer model of the hydraulic transients in the hydropower station is encoded in the FORTRAN programming language. The system as shown in Fig. 4 is divided into 397 reaches. To satisfy Courant's condition for stability of the FD scheme, a time interval of 0.0037s is used. The wave velocities in the pipes are slightly adjusted

X. Yu et al. / Energy 114 (2016) 613e622

so that there is no interpolation error. The characteristic data of the turbine is obtained from characteristic curves as shown in Fig. 5. Two groups of the parameters Td and bt of the governor are selected randomly, where A¼(0.3,8) in the stable region and B¼(0.5,12) in the unstable region. The corresponding eigenvalues are shown in Table 4, which accord with the results in Fig. 6. By using these two groups of governor parameters, the small fluctuation condition is simulated with two different models. As shown in Fig. 7, the variation of dimensionless rotating speed and water level in the surge tank is convergent in both models when the governor parameter A is used. When the governor parameter B is used, the result is divergent in both models, which is shown in Fig. 8. Theoretical analysis is in good agreement with numerical simulations. There are some discrepancies between these two models, because many factors such as water hammer, nonlinear characteristics of the turbine and the effect of the throttle orifice of the surge tank in particular are not considered in the state equations model. If the load change is very small, the results between these two models are almost the same. So it demonstrates that the state equations derived by the proposed method are correct for the GTH system. 5. Conclusions In this work, the state space model of the GTH system under small fluctuation condition is established, and the regular features of the state equations that describe the small fluctuation of the system are investigated. The equations for unsteady flow in the complex pipeline system are conveniently given by using graph theory. By specifying the order of the state variables and using matrix transformation, an innovative method for solving the coefficient matrix of the state equations is established, which is also convenient for programming, and the stable region of the system can be given with the eigenvalue method. The proposed method is used to analyze the stability of the GTH system of an actual hydropower station during small fluctuation, and stability region is given, which reveals the influence of the parameters on system stability. Meanwhile, the small fluctuation in the GTH system is simulated in time domain as well, and theoretical analysis is in good agreement with numerical simulations. This paper provides a novel method for establishing the model of complex GTH system under small fluctuation conveniently, which is the basis for stability analysis and system optimization including the tuning of the governor. Compared with the previous research object, this paper focuses on the whole GTH system of hydropower plants, including the various hydraulic system, mechanical system and electrical system. By using the proposed method, it is convenient to establish the state space model of the complex GTH system and analyze the stability accurately and quickly, which is applicable to consider the effect of the synchronous generator and multi-hydropower plants connection as well. The work carried out in this paper is based on a model using some simplified boundaries, and some practical factors are neglected, such as the effect of the excitation system and electric power system. The authors of this paper intend to continue this work and improve the robustness of the model. Acknowledgments This paper was supported by the National Natural Science Foundation of China (grant numbers: 51379064 and 51409087), the Natural Science Foundation of Jiangsu Province (grant number: BK20130839), and the Fundamental Research Funds for the Central Universities (grant number: 2016B04914).

621

Nomenclature

Acronyms FD finite difference method GTH governor-turbine-hydraulic HVDC high-voltage direct current PID proportional, differential and integral Symbols A, Asu a bp, bt D,D1 e, f, r GD2 H K L M N, NQ n P, PG Q q s1-10, sp

area of pipe, cross-sectional area of surge tank pressure wave celerity temporary speed drop, permanent speed drop equivalent diameter of pipe, turbine diameter coefficient expression moment of inertia of unit head across turbine number of pipeline sections length of pipe number of state variables number of units, number of flow variables rotating speed turbine output, power absorbed by generator discharge dimensionless discharge variation characteristic coefficients of turbine, load self-regulating coefficient T,Td,Tm,Tn torque, dashpot time constant, mechanical starting time, promptitude time constant x dimensionless external load variation Zsu water level of surge tank a head loss coefficient 4 dimensionless rotating speed variation m dimensionless wicket gate opening variation l eigenvalues of coefficient matrix x dimensionless unit's head zin,zout discharge coefficient of throttle orifice t dimensionless wicket gate opening Subscript 11 unit value 0 initial steady state R rated value References [1] Chang XL, Liu XH, Zhou W. Hydropower in China at present and its further development. Energy 2010;35:4400e6. [2] Dorji U, Ghomashchi R. Hydro turbine failure mechanisms: an overview. Eng Fail Anal 2014;44:136e47. [3] Konidaris DN, Tegupoulos JA. Investigation of oscillatory problems of hydraulic generating units equipped with Francis turbine. IEEE Trans Energy Conv 1997;12:419e25. [4] Barelli L, Liucc L, Ottaviano A, Valigi D. Mini-hydro: a design approach in case of torrential rivers. Energy 2013;58:695e706. [5] Yu XD, Zhang J, Hazrati A. Critical superposition instant of surge waves in surge tank with long headrace tunnel. Can J Civ Eng 2011;38(3):331e7. [6] Kishor N, Saini RP, Singh SP. A review on hydropower plant models and control. Renew Sustain Energy Rev 2007;11(5):776e96. [7] Fang HQ, Cheng L, Dlakavu N, Shen ZY. Basic modeling simulation tool for analysis of hydraulic transients in hydroelectric power plants. IEEE Trans Energy Conv 2008;23(3):834e41. [8] Wylie EB, Streeter VL, Suo LS. Fluid transients in systems. Englewood Cliffs, NJ,USA: Prentice-Hall, Inc; 1993. [9] Chaudhry MH. Governing stability of a hydroelectric power plant. Water Power 1970:131e6. [10] Hagihara S, Yokota H, Goda K, Isobe K. Stability of a hydraulic turbinegenerating unit controlled by PID governor. IEEE Trans. Power Apparatus Syst 1979;98:2294e8. [11] Thorne DH, Hill EF. Field testing and simulation of hydraulic turbine governor

622

X. Yu et al. / Energy 114 (2016) 613e622

performance. IEEE Trans Power Appar Syst 1973;92:1183e91. [12] Kishor N, Singh SP, Raghuvanshi AS. Dynamic simulations of hydro turbine and its state estimation based LQ control. Energy Convers Manag 2006;47: 3119e37. [13] Chen DY, Ding C, Do Y, Ma XY, Zhao H, Wang YC. Nonlinear dynamic analysis for a Francis hydro-turbine governing system and its control. J Frankl Inst 2014;351:4596e618. [14] Zhang H, Chen DY, Xu BB, Wang FF. Nonlinear modeling and dynamic analysis of hydro-turbine governing system in the process of load rejection transient. Energy Convers Manag 2015;90:128e37. [15] Zhou JX, Cai FL, Wang Y. New elastic model of pipe flow for stability analysis of the governor-turbine-hydraulic system. J Hydraulic Eng. Asce 2011;137(10):1238e47. [16] Xu BB, Chen DY, Zhang H, Wang FF. Modeling and stability analysis of a fractional-order Francis hydro-turbine governing system. Chaos, Solit Fractals 2015;75:50e61. [17] Guo WC, Yang JD, Yang WJ, Chen JP, Teng Y. Regulation quality for frequency response of turbine regulating system of isolated hydroelectric power plant with surge tank. Electr Power Energy Syst 2015;73:528e38. [18] Yang JD, Wang MJ, Wang C, Guo WC. Linear modeling and regulation quality analysis for hydro-turbine governing system with an open tailrace channel. Energies 2015;8:11702e17. [19] Guo WC, Yang JD, Wang MJ, Lai X. Nonlinear modeling and stability analysis of hydro-turbine governing system with sloping ceiling tailrace tunnel under

load disturbance. Energy Convers Manag 2015;106:127e38. [20] Donaisky E, Oliveira GHC, Santos EAP, Leandro GV, Pena AM, Souza JA. Semiphysical piecewise affine representation for governors in hydropower system generation. Electr Power Syst Res 2016;136:181e8. [21] Martinez-Lucas G, Sarasua JI, Sanchez-Fernandez JA, Wilhelmi JR. Frequency control support of a wind-solar isolated system by a hydropower plant with long tail-race tunnel. Renew Energy 2016;90:362e76. [22] Natarajan K. Robust PID controller design for hydroturbines. IEEE Trans Energy Convers 2005;20(3):661e7. [23] Fang H, Chen L, Shen Z. Application of an improved PSO algorithm to optimal tuning of PID gains for water turbine governor. Energy Convers Manage 2011;52(4):1763e70. [24] Chen ZH, Yuan YB, Yuan XH, Huang YH, Li XS, Li WW. Application of multiobjective controller to optimal tuning of PID gains for a hydraulic turbine regulating system using adaptive grid particle swam optimization. ISA Trans 2015;56:173e81. [25] Chen ZH, Yuan XH, Tian H, Ji B. Improved gravitational search algorithm for parameter identification of water turbine regulation. Energy Convers Manage 2014;78:306e15. [26] Yuan XH, Chen ZH, Yuan YB, Huang YH. Design of fuzzy sliding mode controller for hydraulic turbine regulating system via input state feedback linearization method 2015;93:173e87. [27] Yu XD, Zhang J, Miao D. Innovative closure law for pump-turbines and field test verification. J Hydraulic Eng. Asce 2015;141(3).