Simulation of transient flow in pipeline systems due to load rejection and load acceptance by hydroelectric power plants

ARTICLE IN PRESS International Journal of Mechanical Sciences 52 (2010) 103–115

Contents lists available at ScienceDirect

International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci

Simulation of transient flow in pipeline systems due to load rejection and load acceptance by hydroelectric power plants M.H. Afshar, M. Rohani n, R. Taheri Department of Civil Engineering, Iran University of Science and Technology, Narmak, Tehran, Iran

a r t i c l e in fo

abstract

Article history: Received 29 September 2008 Received in revised form 12 October 2009 Accepted 31 October 2009 Available online 6 November 2009

In this paper, the implicit method of characteristics (IMOC) is extended to simulate the transient flow caused by load variation of hydroelectric power plants. The IMOC was proposed recently by the authors to remove the limitations of the commonly used conventional method of characteristics (MOC) and subsequently used for the simulation of transient flow caused by valve closure and pump failure. In the IMOC, an element-wise definition is used for all devices that may be used in a pipeline system and the corresponding equations are then assembled to form a system of equations that is solved for the unknown nodal heads and flows at each time step. In this paper, the IMOC is extended for the simulation of transient flow in the pipeline system upstream of a hydroelectric power plant. For this, an element-wise definition of the hydroelectric power plants and simple surge tank is proposed and their governing equations are derived in a matrix form required by the IMOC. Then these equations are assembled along with those of other devices of the pipeline system to construct a system of equations that is solved for the unknown head and flow. To further improve the convergence characteristics of the method, an improved formulation of the turbine is also proposed. In the proposed formulation, one of the turbine characteristic parameters, namely turbine net head, is used in addition to turbine head and discharges, to formulate the behavior of the turbine element. The method is applied to a problem of transient flow caused by load rejection and acceptance and the results are presented and compared with those of the explicit MOC. The results show the ability of the proposed method to accurately predict the variations of head and flow in the cases considered. It is also shown that the improved formulation of the turbine is computationally more efficient than the original formulation while producing the same results. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Water hammer Method of characteristics Transient flow Implicit method Power plants

1. Introduction Hydraulic transient events are caused during a change in state, from one steady or equilibrium condition to another. Some of the common events of transient are sudden valve opening or closure, starting or stopping of pumps, mechanical failure of an item, rapid changes in demand condition, etc. The main components of this disturbance are pressure and flow changes that bring about propagation of pressure waves throughout the system. The velocity of this wave may exceed 1000 m/s, which may lead to severe damages. Design and operation of any pipeline system require that the distribution of head and flow in the system is predicted at different operating conditions. Various numerical approaches have been introduced for calculation of the pipeline transients. These include the method

n

Corresponding author. E-mail addresses: [email protected] (M.H. Afshar), [email protected] (M. Rohani), [email protected] (R. Taheri). 0020-7403/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijmecsci.2009.10.014

of characteristics (MOC) [1], wave characteristics method (WCM) [2], finite volume method (FVM) [3], finite element method (FEM), and finite difference method (FDM) [1]. Among these methods, MOC is the commonly used method because of its simplicity and superior performance in comparison with other methods. Koelle and Luvizotto [4] formulated a computer aided operation (CAO) model for the operation of pumped storage plants. The model was based on the use of the MOC to obtain steady, transient, and oscillatory conditions of operation. Kameswara and Eswaran [5] developed a computer program, HYTRAN, using the method of characteristics for the calculation of time dependent head and velocity of fluids in a complex pipeline network for events such as pump failure, load reduction on a turbine, etc. The model considered only a single-phase fluid flow. The analysis had to be modified using a FEM model for the incorporation of structural members such as pipes and penstocks when the effects of fluid–structure interactions (FSIs) were important. Karney and McInnis [6] proposed an explicit algorithm for a general hydraulic element called an external energy dissipator, such as surge tank, relief valve, or reservoir. The method remained explicit even

ARTICLE IN PRESS 104

M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

though friction losses were considered. Ghidaoui et al. [7] provided both a historical perspective and review of the water hammer theory and an overview of recent developments in this field of fluid mechanics. They discussed some of the more complex and fundamental fluid mechanics issues like relation between state equations and wave speeds in single as well as multiphase and multi-component transient flows, various forms of one- and two-dimensional water hammer equations, and governing equations of turbulent water hammer flows by averaging the quasi-two-dimensional plane wave equations. Ramos and Almeida [8] proposed a new technique of analysis based on a dynamic orifice concept. The method sets up the interaction between powerhouse and hydraulic circuits in small hydro schemes by the characteristic parameters of each component. Kochupillai et al. [9] provided a finite element formulation for the fluid–structure interaction (FSI) occurring in pipeline systems mainly due to transient events like rapid valve closure. It was found that there are situations where the structural vibration response increases with the FSIs depending on the relative time scales associated with the structure, flow, and stimulation force. Selek et al. [10] used the method of characteristics to predict the transient behavior of penstock flow of the C - atalan Power Plant in Turkey. Two computational schemes, namely fixed-grid system with space-line interpolation and variable-grid system, were used for the numerical solution. Comparison between the computational results of the transient pressure just upstream of the turbine with those of experimental prototype tests showed reasonable agreement accuracy of these schemes for load rejection and emergency shut down cases. Eker and Grimble [11] simulated the water supply system of Gaziantep city in Turkey, consisting of a sequence of pumping stations that deliver water through pipelines to intermediate storage reservoirs, for management purposes such as flow regulation and cost minimization. The hydraulic models included the nonlinear coupling between flows and reservoir heads. The improved bisection numerical solution method was used to numerically solve the nonlinear equations for the friction coefficient using iterative type of solution method. Katrasnik [12,13] proposed a new algorithm for the simulation of a turbocharger gas turbine and compressor based on extended equations of the method of characteristics. In this method, the inlet and outlet pressure variables as well as gas property variables and concentration can be considered. New equations for the simulation of boundary elements suitable for considering the variation in gas properties were also derived. The improved equations and algorithm guaranteed a much better conservation of mass flux at the expense of slightly increased computational time. Khaleghi et al. [14] presented a onedimensional model for simulation of steady and transient flow through an entire jet engine. A stage-by-stage model of compressor and turbine and simulation of the combustion chamber by a zero dimensional model were also proposed. The method of characteristics was used as the numerical scheme for solving the governing equations. More recently, Afshar and Rohani [15] proposed an implicit method of characteristics (IMOC) as a remedy to the shortcomings of the conventional method of characteristics. In the IMOC, the governing equations of the devices used in a pipeline system were derived in an elementwise manner. These equations were then assembled to form the final system of equations that is solved for the unknown nodal heads and flows. This method was shown to be a versatile method capable of simulating pipeline systems with any arbitrary combination of devices in the system. The performance of the method was demonstrated by some examples of transient flows, including valve closure and pump failure. Application of the IMOC for the simulation of transient flow in pipeline system upstream of the hydroelectric power plants is

considered in this paper. For this, two new devices, namely hydroelectric power plant and surge tank, are introduced and formulated in an element-wise manner and their governing equations are derived in a proper form required by the IMOC. Two basic and improved formulations of the hydropower electric element are also suggested. The proposed methods are used for the prediction of transient flow in a pipeline system consisting of a reservoir, surge tank, and a power plant for two cases of load rejection and load acceptance, and the results are presented and compared with those of the conventional MOC. Experiments indicate the versatility of the IMOC and in particular, the efficiency of the improved power plant formulation compared with the basic one. The explicit method of characteristic is described along with appropriate boundary conditions required for the surge tank and turbine in Section 2. In Section 3, basic concepts of the implicit method of characteristic (IMOC) suggested by Afshar and Rohani [15] are first described for pipe and reservoir elements. Then the underlying concept is extended to drive the governing equations for two new devices, namely surge tank and turbine elements, in the proper form required by the IMOC. An improved formulation of turbine elements is also presented in this section. The efficiency of the improved turbine formulation and application of the IMOC is investigated in Section 4 against a numerical example. And finally the paper is concluded with concluding remarks in Section 5.

2. Conventional explicit MOC The method of characteristics uses the characteristic equations at each internal point of the pipe, shown in Fig. 1, to arrive at the following values of the head and flow at the point under consideration in terms of their values in neighboring point: C þ : QP ¼ CP  Ca HP CP ¼ QA þCa HA  R Dt QA jQA j

ð1Þ

C  : QP ¼ Cn þ Ca HP Cn ¼ QB  Ca HB  RDtQB jQB j

ð2Þ

where R= f/2DA= term of friction factor, Q= discharge, H= pressure head, A= area of the pipe, a= velocity of the pressure wave, f=friction factor of the pipe, Dt =time step size, D= diameter of the pipe, Ca = gA/a, QP =unknown flow at point P at time t + Dt, HP = unknown head at point P at time t+ Dt, QA, QB = flows at neighboring sections of P at the previous time t, and HA, HB = heads at neighboring sections of P at the previous time t. Discretising each pipe of the pipeline system into segments of length Dx, the above equations can be used to calculate the pressure head and flow at all interior points of the pipes. Calculation of flow variables at two end-points of each pipe segment requires appropriate boundary conditions depending on the type of devices used in the system. A pipeline system is a collection of pipes, as the main structure of the system, and various types of devices such as pumps,

t P

t + Δt t

C

+

C

−

B

A Δx

Δx

Fig. 1. Characteristic lines in x–t plan.

x

ARTICLE IN PRESS M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

turbines, check valves, pressure relief valves, surge tanks, etc. In the conventional MOC, the effect of these devices on transient response of the system is considered as boundary conditions for the pipes. In what follows, the boundary conditions imposed by the simple surge tanks and turbines in the pipelines are briefly addressed. A comprehensive discussion of the most common devices can be found in [16]. For a tank located in the middle of a pipeline system, shown in Fig. 2, the following equations can be written: QPi;n þ 1 ¼ QPi þ 1;1 þ Qps

ð3Þ

HPi;n þ 1 ¼ HPi þ 1;1 ¼ zp

ð4Þ

zp ¼ z þ

ð6Þ

Now, other values can be calculated using Eqs. (3)–(5). For a turbine located downstream of a pipeline system, shown in Fig. 3, the appropriate boundary condition can be written as follows [16]: Q2 2gA2t

ð7Þ

where H and Q=pressure head and flow at the entrance to the scroll case, respectively Htail = tail water level above the datum, Hn =net head, and At = cross sectional area of the pressure conduit at the turbine inlet. Calculation of the unknown values of H and Q requires that the net head, Hn, is also calculated. The value of the net head is, in turn, dependent on other characteristics of the turbine such as turbine speed and power through turbine characteristic curves. The characteristic curves of turbine are often graphically defined, as shown in Fig. 4, in terms of the gate opening (t), unit speed (f), Water level at the end of time step

Q ps

Zp Pipe i

2gA

2

Hn

Hp Turbine

Qp H tail

Datum

ð5Þ

Cp  Cn þQs þ ð2As z=DtÞ ¼ Cai þ Cai þ 1 þ ð2As =DtÞ

H ¼ Hn þ Htail 

2

QP

Re servoir

Fig. 3. Turbine at the downstream end of the system [16].

1 Dt ðQps þ Qs Þ 2 As

where QP =discharge at the end of time step, HP = pressure head above the datum, z and zp =height of the liquid surface in the tank above datum at the beginning and end of the time step, respectively, As =horizontal cross sectional area of the surge tank in plane, Qps and Qs =flow through the tank at the end and beginning of the time step, respectively, Dt =time step size, and i represents the pipe. The first equation describes the continuity of flow in the surge tank while the second one considers the pressure head of the pipes connected to the surge tank by neglecting losses at the junction of the surge tank and pipes, and the last equation describes the water surface elevation of the surge tank. Combining the above equations with the positive and negative characteristic line results in the following relation for the head [16]: Hpi;n þ 1

105

Water level at the begining of time step

unit discharge (q), and unit power (p) defined as follows in the SI system:

j¼

DN pffiffiffiffiffiffi 84:59 Hn

q¼

Q pffiffiffiffiffiffi D 2 Hn

i n+1

( i, n + 1 )

QP

i+1,1

P pffiffiffiffiffiffi Hn3

ð8Þ

q ¼ a0 þ a1 j

ð9Þ

where a0 and a1 = constants for the straight line. Substituting Eq. (8) in the above equation leads to [16]: 1=2

a2 Hn

¼ Q  a3

ð10Þ

2

3

where a2 = a0D and a3 = ND a1/84.59. The following relation for the flow is obtained by combining the above equation with the positive characteristic line [16]: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a5  a25  4a4 a6 ð11Þ Q¼ 2a4 in which a4 ¼ Ca =2gA2  Ca =a22 , a5 ¼ 2a3 Ca =a22  1, and a6 ¼   C a2 Cp  Ca Htail  aa 2 3 ; Hn and H are determined from Eqs. (10) 2

and (7), respectively. The next unknown that must be calculated is the rotational speed (N). Variation of the rotational speed (N) is governed by the following ordinary differential equation [16]: Tu ¼ Ttur  Tgen ¼ I

dw 2p dN ¼I 60 dt dt

ð12Þ

where Tu = unbalanced torque, Ttur =turbine torque, Tgen =generator torque, w and N= rotational speed of turbogenerator in rad/s and rmp, respectively, and I=total moment of turbine and generator. Integration of this equation on a time step followed by some mathematical manipulations leads to the following equation, which can be used to calculate the rotational speed at each time step of the computation [16]: ( ðN n Þ2 þ 182; 380

Dt I

Pipe i + 1 QP ,

D2

where D=diameter of runner, N=rotational speed, and P =power output. A linear extrapolation from the q  f curve is usually used to arrive at the following relations:

Nn þ 1 ¼ Z

p¼

" n nþ1 þ Ptur Þ 0:5ðPtur

0:5

Zgen

#)0:5 n nþ1 ðPgen þ Pgen Þ

ð13Þ where Pgen = generator load, Ptur =power developed by turbine, and

Zgen = generator efficiency.

( i + 1,1 ) Datum

Fig. 2. Surge tank in the middle of the system.

Variation of the gate opening during transients is governed by the following four differential equations in terms of four variables of actuator rate (va), dashpot saturation (et), gate

ARTICLE IN PRESS 106

M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

Gate opening ,τ (%)

100 Gate opening ,τ (%)

0

80

0.05 100 90

80

0

φ Unit speed,

90

Unit power, p

Unit discharge, q

0.5

φ Unit speed,

Fig. 4. Typical characteristic curves of turbine [16].

servomotor rate (vd), and gate opening (t) [16]:   dva 1 det 1 dv ¼ ¼ ðnref  n  et  sva Þ; DTr a  et dt dt dt Ta Tr  dt dvd 1 ¼ ¼ ks vd k ðva  tÞ  vd ; dt dt Td d

ð14Þ

where n =N/Nr, nref =1+ st0, Ta = actuator time constant, Tr = temporary time constant, d =temporary speed droop, s = permanent speed droop, Td =distributing valve time constant, kd =distributing valve gain, ks = gate-servomotor gain, Nr = synchronous turbogenerator speed, and t0 = initial steady state gate opening. As is clear from the above equations, the four parameters turbine head, flow, rotational speed, and gate opening are unknown. For the calculation of these unknown values at each time step, the following procedure is suggested by [16], which is schematically shown in Fig. 5: 1- estimate the values of N; t, and Hn at the end of ith time step from the known values of the previous three time steps, i-1, i-2, and i-3, by parabolic extrapolation; 2- calculate the turbine characteristic parameters (f,q) from estimated values of the previous step from Eq. (8); 3- calculate Q from Eqs. (11) and (10); 4- the estimated value of j is calculated using the estimated value of N and computed value of Hn; 5- calculate Ptur from the turbine characteristic data and N from Eq. (13); 6- if N  N 40:001Nr , N is assumed equal to N and the above procedure is repeated; 7- use Eq. (14) for calculation of t; 8- if jt  tj 4 0:005, t is assumed equal to t and the above procedure is repeated. To prevent unlimited repetition of iterations, a counter is often introduced in both iterative loops [16].

3. Implicit method of characteristics Different methods have been proposed and successfully used for transient flow simulation of pipeline systems. Of these, the well-known method of characteristics (MOC) has been mostly used for the transient flow analysis and in particular for transient flow simulation in hydropower plants [16,17]. This method has many attractive features; simplicity is the most important one. The method uses governing equations of the transient flow on characteristic lines to arrive at the algebraic equations to be solved for the head and flow at each interior point of the pipes. The devices used in the system are taken into account by

considering their governing equations as boundary conditions for the end points of the pipes. Computer implementation of the conventional explicit MOC, however, faces some serious limitations, as application of this method requires that only one device is located between any two pipes. Otherwise, a new boundary condition has to be derived for each combination of devices that may be used in the pipeline system. Furthermore the boundary conditions derived for different devices are very much dependent on whether these devices are located at the end of the system or within the system. To further complicate the situation, these boundary conditions change depending on whether these devices are located at the upstream or downstream end of the system. An implicit method of characteristics (IMOC) was proposed by Afshar and Rohani [15] to remove the shortcomings of the conventional MOC. In this method, an element-wise definition was used for all devices of a pipeline system and the corresponding equations were derived in a proper form to be used in the IMOC. These equations were independent of the location of devices and flow direction. Then they were assembled to form the system of equations to be solved for the unknown head and flow at nodes. The formulation was shown to accommodate any arbitrary combination of the devices in a pipeline system. This method was successfully used for the simulation of transient flows caused by sudden valve closure and pump failure [15] in the pipeline systems. Here, this method is extended to consider devices such as surge tanks and turbines and is subsequently used for the simulation of transient flow in the pipeline containing hydroelectric power plant. Any hydropower electric plant system consists of at least a reservoir with fixed head upstream, a hydroelectric power plant downstream, a set of pipes conveying water from upstream reservoir to the power plant, and usually a surge tank located upstream of the power plant to damp the water hammer that may occur due to the action of the power plant. In what follows, two main components of a hydropower electric plant system, namely surge tank and turbine, are formulated in the proper form to be used in the IMOC analysis of the pipeline system. Before that, however, the formulations of other components such as pipe and reservoir as suggested in [15] are briefly reviewed for the sake of completeness. More detailed discussion on the development of pipe and reservoir elements along with other type of devices can be found in [15].

3.1. Pipe Consider a pipe segment defined by two end points i and j shown in Fig. 6. The proposed formulation for pipes starts with characteristic equations. Here, a second order central differencing

ARTICLE IN PRESS M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

107

Call From Main Program Return Parabolically Extrapolate N , τ And Comput e H n From Eqs. 1 And 3

Print τ , N , H n , Q & P Yes

Compute Q From N And H n No

Step For Print ing Values

Return Det ermine Q From Eq. 7 And H n From Eq. 6 τ =τ

Stop

Comput eφ From N And H n

N =N Print " Iterations For τ Failed"

Det ermine Ptur From T urbine Characteristics Dat a

No N =N No. Of Iterat ions < 30?

Compute N FromEq. 9

No

No. Of Iterat ions > 30?

Yes

No

N−N

Yes No

τ − τ p 0.005

p 0.001

Nr Yes

Print "Iterat ions For N Failed"

Compute τ From Governor Subrout ine

Stop

Fig. 5. Flowchart of boundary condition for a Governed Francis Turbine [16].



t

( n + 1 ) t + Δt

C

Δt t+ 2

(n) t

−

C

+

i

An analogous treatment allows for writing the equivalent equations for the negative characteristic:

j

x

Fig. 6. Pipe element used in IMOC.

C :

is used to arrive at the required algebraic form of the positive characteristic: Cþ :

dQ dH þ Ca þRQ jQ j ¼ 0 ðn þ ð1=2ÞÞ dt dt

 1 1 R Dt Qjn þ 1 signðQjn þ 1 þQin Þ þ R Dt Qin signðQjn þ 1 þ Qin Þ Qjn 4 2 1 þ 1 þCa Hjn þ 1 ¼ Qin þ Ca Hin  R Dt ðQin Þ2  signðQjn þ 1 þQin Þ 4 ð16Þ

1þ

ð15Þ

dQ dH  Ca þRQ jQ j ¼ 0 ðn þ ð1=2ÞÞ dt dt

ð17Þ

 1 1 1 þ R Dt Qin þ 1 signðQin þ 1 þQjn Þ þ R Dt Qjn SignðQin þ 1 þ Qjn Þ Qin 4 2 1 þ 1  Ca Hin þ 1 ¼ Qjn  Ca Hjn  R Dt ðQjn Þ2  signðQin þ 1 þ Qjn Þ 4 ð18Þ

ARTICLE IN PRESS 108

M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

where Hi, Hj = values of the pressure head at the end points of pipe segment, Qi, Qj = values of the flow at the end points of pipe segment, and n represents the time step. Eqs. (16) and (18) totally describe the hydraulic behavior of a pipe segment in terms of the four unknown parameters Qin þ 1 , Qjn þ 1 , Hin þ 1 , and Hjn þ 1 of the pipe element. These equations, however, are clearly non-linear and their solutions require a linearization scheme. Here a Newton–Raphson formulation is used to lead the following linear system of equations: e

S epipe DX epipe ¼ bpipe

with

DX Tpipe ¼ ½DQi ; DQj ; DHi ; DHj n þ 1

3.2. Surge tank element The equations governing the behavior of a surge tank, shown as a hydraulic element in Fig. 7, can be defined as follows: Qin þ 1 ¼ Qjn þ 1 þ Qsn þ 1

Hin þ 1 ¼ zn þ 1 þ ks

Qsn þ 1 Qsn þ 1 2 2gAs

ð23Þ

Hjn þ 1 ¼ zn þ 1 þ ks

Qsn þ 1 Qsn þ 1 2gA2s

ð24Þ

ð19Þ

where DQ=Qm + 1  Qm, DH=Hm + 1  Hm, and m is the nonlinear e iteration index; S epipe and bpipe are the element matrices defined as " # 0 C 0 Ca S epipe ¼ D 0 Ca 0

ð22Þ

dz Qs ¼ As dt

ð25Þ

2

3 1 1 Qjm þ R Dt ððQin Þ2 þ ðQjm Þ2 Þ signðQjm þ Qin Þ þ R Dt Qjm Qin signðQjm þQin Þ þ Ca Hjm  Qin  Ca Hin 7 6 4 2 e 7 bpipe ¼  6 4 m 1 5 1 Qi þ R Dt ððQim Þ2 þ ðQjn Þ2 Þ signðQim þ Qjn Þ þ R Dt Qjn Qim signðQim þQjn Þ  Ca Him  Qjn þ Ca Hjn 4 2

with C ¼ 1þ

1 1 R Dt Qjm signðQjm þ Qin Þ þ R Dt Qin signðQjm þ Qin Þ 2 2

and D ¼ 1þ

1 1 R Dt Qim signðQim þ Qjn Þ þ R Dt Qjn signðQim þQjn Þ 2 2

ð20Þ

The system of equations describing the hydraulic behavior of other devices such as reservoir, valve, junction, and pump can also be derived using the same idea described for the pipe element [15]. Here only the system of equations governing the hydraulic behavior of the reservoir element is reported from [15]. More detailed discussion of this element can be found in [15]. Reservoir element: The stiffness matrix and right hand side vector of a reservoir element were found to be [15]: e

S eres DX eres ¼ bres

with

DX Tres ¼ ½DQi ; DQj ; DHi ; DHj n þ 1

2

S eres

1 6 ¼ 4 2ð1 þ kres Þ ðQ m þ Q m ÞsignðQ m þ Q m Þ i j i j 8gA2res

2 6 e bres ¼  4

Him



Hjm

m+1

1

0

2ð1 þ kres Þ m ðQi þQjm ÞsignðQim þ Qjm Þ 8gA2res

1

3

Qim  Qjm

7 ð1 þkres Þ m 5  ðQi þQjm Þ2 signðQim þ Qjm Þ 2 8gAres m+1

m

where Qs =flow through the tank, ks = coefficient of entrance loss of the surge tank, z =height of the liquid surface in the tank above datum, As = horizontal cross sectional area of the surge tank in plane, Hi and Hj = values of the pressure head at the first and second node of the surge tank, Qi and Qj = values of the flow at the first and second node of the surge tank, respectively, and n represents the time step. The first equation describes the continuity of flow in the surge tank while the second and third equations define the surge tank element nodal pressure heads and the last equation describes the water surface elevation in the surge tank. Now Eqs. (22) and (24) are used to derive the required element matrices. These equations, however, cannot be used in their present form to derive the surge tank matrices since the elevation of the liquid surface in the tank (zn + 1) and the flow through the tank (Qsn þ 1 ) are not known. For this, a second order central differencing is first used on Eq. (25) to arrive at the following equation for the water surface elevation in the tank:

zn þ 1 ¼ zn þ ð21Þ

where DQ=Q  Qm, DH= H  H , m is the nonlinear iteration index, kres = coefficient of entrance loss of the reservoir, Ares =area of the reservoir orifice taken as the average cross sectional area of the pipes connected to the reservoir, Hi and Hj = values of the pressure head at the first and second node of the reservoir, Qi and Qj = values of the flow at the first and second node of the reservoir, respectively, and n represents the time step. The idea behind the development of the IMOC is now used in this paper to formulate the hydraulic behavior of two other type of devices, namely surge tank and turbine, which are commonly used in the hydroelectric power plant.

0

3

7 1 5

1 Dt n ðQ þQsn þ 1 Þ 2 As s

ð26Þ

where Dt is the duration of time step. Replacing zn + 1 from Eqs. (26) and (22) into Eqs. (23) and (24) leads to the required form of the governing equations for the

Qs i

j

z

Datum Fig. 7. Surge tank element used in IMOC.

ARTICLE IN PRESS M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

surge tank: Hin þ 1 ¼ zn þ

Hin þ 1  Hjn þ 1 ¼ Hnn þ 1 

Dt

ðQin þ 1  Qjn þ 1 þQin  Qjn Þ

2As ks þ ðQin þ 1  Qjn þ 1 Þ Qin þ 1  Qjn þ 1 Hjn þ 1 2 2gAs Dt n þ 1 ¼ zn þ ðQ  Qjn þ 1 þ Qin  Qjn Þ 2As i ks þ ðQ n þ 1  Qjn þ 1 Þ Qin þ 1  Qjn þ 1 2gA2s i

ð27Þ

These equations can now be written in a matrix form as defined by Eq. (19) using a Newton–Raphson method of linearization to get the following stiffness and right hand side matrices for the surge tank element: e

S etan k DX etan k ¼ btan k ;

DX Ttan k ½DQi ; DQj ; DHi ; DHj n þ 1

ð28Þ

where DQ= Qm + 1 Qm, D =H= Hm + 1 Hm, m is the nonlinear iteration index, and:

 B B 1 0 S etan k ¼ B B 0 1

109

ðQjn þ 1 Þ2

ð33Þ

2gA2t

These equations, however, are not in the proper final form since the net head is not known. The value of the net head is a function of dimensionless parameters of the turbine through turbine characteristic curves as explained in previous section. These dependencies are usually defined by two set of equations graphically, shown in Fig. 4. Here, the following equations of the general forms are used to explain the formulation for turbine simulation: Fðt; j; qÞ ¼ 0;

Gðt; j; pÞ ¼ 0

ð34Þ

These equations coupled with Eqs. (13) and (14), describing the variation of rotational speed and the gate opening at each time step, respectively, are used to derive the turbine element governing equations in a proper form to be used in IMOC. Using the definition of net head in terms of the unit discharge (q) from Eq. (8) into Eqs. (32) and (33) leads to the following relations: Hin þ 1  Hjn þ 1 ¼

ðQin þ 1 Þ2 ðqD2 Þ2



ðQin þ 1 Þ2 2gA2t

ð35Þ

3 Dt m Dt n ks m n m m 2 m m m n  z  ðQ  Q Þ  ðQ  Q Þ  ðQ  Q Þ signðQ  Q Þ H i j i j i j i j 7 6 i 2As 2As 2gA2s 7 6 e btan k ¼  6 7 D t D t k 5 4 m s n n m m 2 m m n Hj  z  ðQi  Qj Þ  ðQi  Qj Þ  ðQi  Qj Þ signðQi  Qj Þ 2 2As 2As 2gAs 2

with

Dt

2ks B¼ þ ðQ m  Qjm ÞsignðQim  Qjm Þ 2As 2gA2s i

ð29Þ

Hin þ 1  Hjn þ 1 ¼

ðQjn þ 1 Þ2 ðqD2 Þ2

ðQjn þ 1 Þ2

ð36Þ

2gA2t

A Newton–Raphson formulation can now be used to arrive at the following stiffness and right hand side matrices for the turbine element:

3.3. Basic turbine element

e

For a turbine element, shown in Fig. 8, two equations of turbine head and continuity of flow are used as the starting point to define the behavior of the turbine element in terms of nodal head and flow parameters as follows: Hin þ 1  Hjn þ 1 ¼ Hnn þ 1 



ðQin þ 1 Þ2 2gA2t

ð30Þ

Qin þ 1 ¼ Qjn þ 1

S eturb DX eturb ¼ bturb ; m+1

where DQ=Q index, and

DX Tturb ¼ ½DQi ; DQj ; DHi ; DHj n þ 1 m

 Q , DH=H

m+1

ð37Þ

m

 H , m is the nonlinear iteration

3 ðQ m Þ2 ðQ m Þ2 6 Him  Hjm  i 2 þ i 2 7 6 2gAt 7 ðqD2 Þ 7 6 e bturb ¼  6 7 6 ðQjm Þ2 ðQjm Þ2 7 5 4 Hm  Hm  þ i j 2gA2t ðqD2 Þ2 2

S eturb ¼



C

0

1

1

0

D

1

1

 ;

ð31Þ

where Hn = the net head, At =cross sectional area of the pressure conduit at the turbine inlet, Hi and Hj = values of the pressure head at the first and second node of the turbine, and Qi and Qj = values of the flow at the first and second node of the turbine, respectively. These equations are first combined to lead to the following proper form defined as Hin þ 1  Hjn þ 1 ¼ Hnn þ 1 

ðQin þ 1 Þ2 2gA2t

ð32Þ

Hydraulic Grade Line

Hi

Hj i

j

Datum Fig. 8. Turbine element used in IMOC.

C ¼ 2Qim

1 ðqD2 Þ2

þ

! 1 ; 2gA2t

D ¼ 2Qjm

1 ðqD2 Þ2

þ

1 2gA2t

! ð38Þ

Calculation of these matrices at each time step of the calculation, however, requires that the value of unit discharge (q) is known. For this, the following procedure is suggested and used for the calculation of unit discharge (q) at each nonlinear iteration of the solution procedure: 1. Use Eq. (30) to calculate an initial estimate of the net head Hn using the most recent values of nodal head and flows. At the start of the nonlinear iterations, the values of nodal head and flows at the previous time step are used to start the computation. 2. Use the power characteristic curve G(t,f,p) to calculate the unit power (p) using the most recent values of N and t. The calculated unit power is then used to update the value of turbine power Ptur. 3. Eq. (13) is used to update the rotational speed (N). 4. The set of ordinary differential equations 14 is solved using the known value of rotational speed (N) to update the value of gate

ARTICLE IN PRESS 110

M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

opening (t). These equations are solved here using the fourthorder Runge–Kutta method. 5. Use the discharge characteristic curve F(t,f,q) to calculate the unit discharge (q) using the most recent values of N and t. 6. Definition of unit discharge is used to update the value of net head (Hn). 7. Go to step 2 if this process is not converged. Experiments show that even one iteration of the above procedure is enough for the convergence of the whole system of nonlinear equations. The convergence of the system of nonlinear equations, however, is affected by the accuracy of the above calculation. More nonlinear iterations are required at each time step if the above procedure is not converged to the required accuracy. 3.4. Improved turbine element A second formulation of the turbine is also possible for improved performance of the IMOC. In this formulation, the stiffness equation of the turbine element is derived in terms of five parameters, i.e.; four parameters of head and flow at two end points and the net head of the turbine. For this, the definition of unit discharge as given in Eq. (8) is also included in the set of turbine governing equations to lead to the following equations: Hnn þ 1 ¼

ðQin þ 1 Þ2

ð39Þ

ðqD2 Þ2

Hin þ 1  Hjn þ 1 ¼ Hnn þ 1 

Hin þ 1  Hjn þ 1 ¼ Hnn þ 1 

ðQin þ 1 Þ2 2gA2t

ð40Þ

ðQjn þ 1 Þ2

ð41Þ

2gA2t

Eqs. (39)–(41) are employed to arrive at the following stiffness and right hand side matrices for the turbine element using a Newton–Raphson linearization method: e

e

S eturb Dturb ¼ bturb ; m+1

DTturb ¼ ½DQi ; DQj ; DHn ; DHi ; DHj n þ 1 m

m+1

m

where DQ=Q  Q , DH=H  H , DHn ¼ nonlinear iteration index, and 3 2 2Qim 0 1 0 0 7 6 7 6 ðqD2 Þ2 7 6 7 6 2Qim 7 6 e 0 1 1 1 S turb ¼ 6 7 2 7 6 2gAt 7 6 m 7 6 2Qj 5 4 0 1 1 1 2 2gAt 3 2 ðQim Þ2 m 7 6 Hn  7 6 ðqD2 Þ2 7 6 6 m 2 7 ðQ Þ 7 6 i m m m e H  H  H þ 7 6 n bturb ¼  6 i j 2 7 2gA t 7 6 7 6 6 ðQjm Þ2 7 5 4 m Hi  Hjm  Hnm þ 2gA2t

Hnm þ 1



ð42Þ Hnm ,

m is the

procedure used in the conventional MOC defined in the previous section. This is clearly evident from the procedure outlined in the previous section for the calculation of turbine flow and net head, which requires convergence check on two internal iterations which are very much prone to divergence. The efficiency and effectiveness of the proposed methods are tested in the next section.

4. Numerical application The efficiency of the proposed models is examined here with two numerical examples of load acceptance and load rejection. The example is the case of a transient flow in a power plant system shown in Fig. 9 solved by the proposed IMOC and the proposed turbine element. Two cases of load acceptance and load rejection of the generator are considered to be the causes of transient. These examples were proposed and investigated by Chaudhry [16] using a conventional MOC. The power plant system consists of two reservoirs, two pipes, a surge tank, and a turbine with the following characteristics: The length of the first pipe is 1964 m, diameter of the first pipe is 5.442 m, friction factor of the first pipe is 0.01144, and the wave velocity in the first pipe is considered to be 146 m/s. The second pipe has a length equal to 146 m, a diameter of 5.442 m, a friction factor of 0.01144, and a wave velocity of 146 m/s. The horizontal cross sectional area of surge tank is assumed to be 148.8 m2. The following values are used for the governing constants: Td =0.05, kd = 10, ks =0.2, Ta = 0.1, and s =0.03 [16]. The characteristic data of the turbine of the Francis type can be found in [16]. The load acceptance is assumed to occur from 10 to 25 MW for 10 s and the load rejection, which causes more severe pressure surges, is considered from 25 to 10 Mw in 10 s. Fig. 10 shows the computational layout of the power plant system in which all the pipes and devices are represented in an element-wise manner required by the proposed IMOC. Figs. 11–13 show the variation of the net head, discharge, and rotational speed of turbine during 3000 seconds obtained by IMOC in two cases. Figs. 14 and 15 also show the variation of the gate opening and water elevation in the surge tank located upstream of the turbine for 3000 s obtained by IMOC in the cases of load acceptance and rejection. This problem is also solved using the explicit MOC, so that a comparison can be made with the results obtained using the implicit method. Fig. 16 shows the relative difference between the turbine’s net head obtained with the implicit and the explicit MOC for load acceptance. This figure indicates a good agreement between the results of the explicit and the implicit MOC. Note that the maximum relative difference between the two results for load acceptance is 0.00136, mounting to a maximum difference of 0.136 m. The authors, however,

Surge tank

ð43Þ

Re servoir

Pipe 1

Pipe 2

Turbine

Re servoir

Water Level in surge tan k Fig. 9. Power plant system of the example.

The same procedure outlined before is used for the calculation of unit discharge (q) at each nonlinear iteration. The modified turbine formulation is expected to lead to faster convergence of the method compared with the basic formulation. The method, however, leads to an increase in the number of system unknowns, one for each turbine, which is negligible considering the improved convergence characteristics of the resulting model. The Proposed formulations can be seen to be theoretically more robust than the

Pipe 1

Pipe 2

Turbine Re servoir

Surge Tank

Re servoir

Fig. 10. Computational layout of the power plant system of the example defined as an assembly of appropriate elements used in IMOC.

ARTICLE IN PRESS M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

100

111

Load Acceptance Load Rejection

Hn (m)

98

96

94

92

90 0

500

1000

1500

2000

2500

3000

Time (s) Fig. 11. Net head vs time obtained by IMOC for load acceptance and rejection.

40

30

3

Flow (m /s)

35

25 20 15

Load Acceptance Load Rejection

10 0

500

1000

1500

2000

2500

3000

Time (s) Fig. 12. Turbine flow vs time obtained by IMOC for load acceptance and rejection.

140

Load Acceptance Load Rejection

N (rpm)

135

130

125

120

115 0

500

1000

1500

2000

2500

Time (s) Fig. 13. Rotational speed of turbine vs time obtained by IMOC for load acceptance and rejection.

3000

ARTICLE IN PRESS 112

M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

0.6

Load Acceptance

Gate Opening (%)

Load Rejection 0.5

0.4

0.3

0.2 0

500

1000

1500

2000

2500

3000

Time (s) Fig. 14. Gate opening vs time obtained by IMOC for load acceptance and rejection.

150

Load Acceptance

Tank Water Level (m)

Load Rejection 148

146 144

142

140 0

500

1000

1500

2000

2500

3000

Time (s) Fig. 15. Elevation of water in the surge tank vs time obtained by IMOC for load acceptance and rejection.

0.14 Comparison of the IMOC and explicit MOC

Relative Difference (%)

0.12

Comparison of the basic and improved formulation

0.1 0.08 0.06 0.04 0.02 0 0

1

2

3

4

5

6

7

8

9

10

Time (s) Fig. 16. Relative difference between turbine net heads obtained by the IMOC and explicit MOC and basic and improved turbine formulation for load acceptance.

ARTICLE IN PRESS M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

0.12

Comparison of the IMOC and explicit MOC Comparison of the basic and improved turbine formulation

0.1

Relative Difference (%)

113

0.08

0.06

0.04

0.02

0 0

2

4

6

8

10

Time (s) Fig. 17. Relative difference between turbine net heads obtained by the IMOC and explicit MOC and basic and improved turbine formulation for load rejection.

12

Number of Iterations

10 8 6 4 2

Basic Turbine Formulation (Rejection) Improved Turbine Formulation (Rejection-Acceptance) Basic Turbine Formulation (Acceptance)

0 0

20

40

60

80

100

Time step Fig. 18. Number of nonlinear iterations required by the improved and basic turbine formulation at each time step for load acceptance and rejection.

1000

Basic Turbine Formulation (Rejection)

Total Number of Iterations

900

Improved Turbine Formulation (Rejection-Acceptance)

800

Basic Turbine Formulation (Acceptance)

700 600 500 400 300 200 100 0 0

20

40

60

80

100

Time step Fig. 19. Total number of nonlinear iterations required by the improved and basic turbine formulation for load acceptance and rejection.

ARTICLE IN PRESS 114

M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

Total Required Time (milli second)

1000

Basic Turbine Formulation (Rejection) Improved Turbine Formulation (Rejection-Acceptance) Basic Turbine Formulation (Acceptance)

800

600

400

200

0 0

10

20

30

40

50

60

70

80

90

100

Time step Fig. 20. Total CPU time required by the improved and basic turbine formulation for load acceptance and rejection.

believe that the solutions obtained by the proposed methods are more accurate within the limitations of the turbine formulation used. This problem is again solved using the improved turbine formulation described earlier. The relative difference between the turbine’s net head obtained with the improved turbine element and those of the basic turbine element is also shown in Fig. 16, indicating that the solutions obtained by both turbine elements are virtually the same. Note that the maximum absolute difference in the solutions is just about 0.03 m. Fig. 17 shows the same results as those of Fig. 16 for load rejection. This figure shows the relative difference of the turbine’s net head obtained with the implicit and the explicit MOC, indicating a maximum relative difference of 0.001 amounting to a maximum difference of 0.1 m between these results. These results clearly show the ability of the proposed method to accurately predict the pipeline response to load rejection. This load rejection problem is also solved using the improved turbine formulation described earlier so that a comparison can be made with the results of the basic turbine formulation. Fig. 17 also shows the relative difference of the turbine’s net head obtained using the basic and improved turbine formulations. This comparison shows that the maximum absolute difference in the solutions is just about 0.0133 m indicating that both solutions are virtually the same. The improved formulation, however, requires much less computational effort represented by the fewer number of iterations required to solve the nonlinear equations of the pipeline system at each time step. This is shown in Fig. 18, where it is seen that the number of nonlinear iterations required by the basic formulation gradually increases during load acceptance and load rejection and stays nearly constant once total load is accepted or rejected. The improved formulation, however, requires the same number of nonlinear iterations, nearly one third of the first one, at all times of the simulation process for both cases of load rejection and acceptance. These results clearly show the efficiency of the improved turbine element compared with the basic one. Fig. 19 compares the variation of the total number of nonlinear iterations required by the proposed improved turbine formulation to those of the basic formulation for two cases of load acceptance and rejection. It is observed that the computational effort required by the improved turbine formulation represented by the total number of nonlinear iterations is about 33 percent of that required by the original formulation for load acceptance, and about 42 percent for load rejection. And finally, Fig. 20 compares the total CPU time on a 2 MHz Pentium 4 required by the original and improved turbine formulation for two cases of load

acceptance and rejection at different time steps of the transient calculation. This Figure shows that the improved turbine formulation on average requires nearly 55 percent of the CPU time required by the basic formulation, indicating the higher efficiency of the method.

5. Conclusion An implicit method of characteristics is used in this paper to simulate the response of a pipeline system upstream of power plants to two cases of load acceptance and load rejection. This method is recently proposed by the authors as a remedy to the shortcomings and limitations of the conventional method of characteristics (MOC). An element-wise definition is proposed for all the devices commonly used in the system and in particular, surge tank and turbine, and the corresponding equations are derived in an element-wise manner. The equations are then assembled to form the system of equations to be solved for the unknown nodal heads and flows. An improved turbine formulation is also proposed to improve the efficiency of the proposed IMOC. In the improved turbine formulation, the stiffness equation of the turbine element is derived in terms of five parameters, i.e. four parameters of head and flow at two end points and one of the characteristic parameters of the turbine, net head, leading to faster convergence of the nonlinear iterations at each time step. The efficiency and effectiveness of the IMOC and proposed turbine formulations are tested against a numerical example of load acceptance and rejection and the results are presented and compared with those of the explicit MOC. The results show the ability of the IMOC and proposed surge tank and turbine elements to accurately predict the system transient response to load acceptance and rejection and also indicate improved efficiency of the proposed turbine formulation.

References [1] Izquierdo J, Iglesias PL. Mathematical modeling of hydraulic transients in simple systems. Mathematical and Computational Modeling 2002;35(7): 801–812. [2] Wood DJ, Lindireddy S, Boulos PF, Karney B, Mcpherson DL. Numerical methods for modeling transient flow. Journal of American Water Works Association (AWWA) 2005;97(7):104–15. [3] Zhao M, Ghidaoui M. Godunov-type solutions for water hammer flows. Journal of Hydraulic Engineering, ASCE 2004;130(4):341–8.

ARTICLE IN PRESS M.H. Afshar et al. / International Journal of Mechanical Sciences 52 (2010) 103–115

[4] Koelle E, Luvizotto E. Operational simulation of pumped storage plant, including transient and oscillatory phenomena. Computer methods in water resource II 1991:317–28. [5] Kameswara CV, Eswaran K. Pressure transients in incompressible fluid pipeline networks. Nuclear Engineering and Design 1999;188(1):1–11. [6] Karney BW, McInnis A. Efficient calculation of transient flow in simple pipe network. Journal of Hydraulic Engineering 1992;118(7):1014–30. [7] Ghidaoui MS, Zhao M, McInnis DA, Axworthy DH. A review of water hammer theory and practice. Applied Mechanics Reviews 2005;58:49–76. [8] Ramos H, Almeida AB. Parametric analysis of water hammer effects in small hydro power. Journal of Hydraulic Engineering, ASCE 2002;128(7):689–96. [9] Kochupillai J, Ganesan N, Padmanabham C. A new finite element formulation based on the velocity of the flow for water hammer problem. International Journal of Pressure Vessels and Piping 2005;82(1):1–14. ¨ MS, Selek Z. Comparison of computed water hammer [10] Selek B, Kirkgoz pressures with test results for the C - atalan power plant in Turkey. NRC Research Press Web Site 2004:78–85.

115

[11] Eker I, Grimble MJ. Operation and simulation of city of Gaziantep water supply system in Turkey. Journal of Renewable Energy 2003;28:901–16. [12] Katrasnik T. A novel algorithm for the simulation of an automotive turbocharger turbine. Journal of Mechanical Engineering 2006;52(1):15–25. [13] Katrasnik T. Improved model to determine turbine and compressor boundary conditions with the method of characteristics. Journal of Mechanical Science 2006;48(5):504–16. [14] Khaleghi H, Boroomand M, Tousi AM. A stage by stage simulation and performance prediction of gas turbine. Proceedings of the 43rd AIAA aerospace sciences meeting and exhibit 2005:59–68. [15] Afshar MH, Rohani M. Water hammer simulation by implicit method of characteristic. International Journal of Pressure Vessels and Piping 2008;85:851–9. [16] Chaudhry HM. In: Applied hydraulic transients. New York: Van Nostrand Reinhold; 1987. [17] Swaffield JA, Boldy A. In: Pressure surge in pipe and duct systems. UK: Gower Press; 1993.