ISSN 00406015, Thermal Engineering, 2011, Vol. 58, No. 13, pp. 1094–1098. © Pleiades Publishing, Inc., 2011. Original Russian Text © E.A. Tairov, A.A. Levin, V.F. Chistyakov, 2011, published in Izvestiya RAN. Energetika.

Using the Theory of Hydraulic Circuits in Simulating Thermal Power Installations E. A. Tairov, A. A. Levin, and V. F. Chistyakov Received May 17, 2010

Abstract—Specific features relating to the motion of media in the paths of thermal power stations are con sidered from the viewpoint of the hydraulic circuit theory. It is shown that difficulties arise when an attempt is made to apply traditional approaches for calculating flow distribution using the methods of the hydraulic circuit theory. A calculation scheme is proposed the use of which makes it possible to reduce the initial prob lem statement to a classic scheme according to the method of nodal pressures. Test results are presented that demonstrate high reliability of the developed method for calculating flow distribution. DOI: 10.1134/S004060151113009X

INTRODUCTION Successful application of decomposition principles in modeling complex power installations, separated description of slow heat and masstransfer processes and fast hydrodynamic processes are factors that prompt specialists to develop and use versatile proce dures for constructing mathematical models. It is of interest to develop a generalized algorithm for calcu lating flow distribution in the process paths of a power installation using the methods of the hydraulic circuit (HC) theory [1]. In [2], it is pointed out that the math ematical description of elements composing the ther mal circuits of thermal power stations (TPSs) lacks uniqueness. The authors draw a conclusion that there is a need to unify the mathematical description of indi vidual components comprising the process diagram of a TPS for making it possible to develop universal com putation codes. Hence, there is a need to develop a generalized algorithm for calculating flow distribution in the steam–water path of a power installation, in which water and steam flow through channels that have different cross sections and lengths and contain

flow connection and separation nodes, stop and con trol members, and sources of head (pumps). The pro cess diagram of a power unit contains hundreds of nodes and branches, which generates the need to solve systems of equations with a large dimension with respect to flowrates and pressures, and the methods of the HC theory are aimed at solving such systems of equations. THE STANDARD STATEMENT OF A FLOW DISTRIBUTION PROBLEM Flow distribution calculations are usually carried out under the assumption that friction has a predomi nant influence on pressure drop, due to which it becomes possible to use essentially simplified mathe matical models for describing the motion of medium in the lengthy process paths of power installations [1, 3]. In this case, an HC can be represented as a graph consisting of m nodes and n branches and described using the matrix of coupling A:

⎧0, if the branch j has no connection with the node i; ⎪ Aij = ⎨1, if flow through the branch j exits from the node i; ⎪⎩− 1, if flow through the branch j enters into the node i;

where i = 1, m, j = 1, n. We introduce the following quantities: the vector of flowrates through the branches D = (d1, d2,..., d m), vec tor of pressures at the nodes P = ( p1, p2,..., pm ), vector of pressure differences across the branches Y = (y1, y2,..., y m ), vector of flowrates at the nodes Q = (q1, q2,..., qm ) , vector of the hydraulic resistances of the

branches S = (s1, s2,...s n), and vector of effective heads in the branches H = (h1, h2,..., h n) . The components of the vector of flowrates are interconnected by the rela tion ∑ q j = 0. The flowrates through the HC and the pressures in it are interconnected by the following equations [1, 3]:

1094

T

AD = Q; A P = Y ; Y + H ≡ SDD ,

(1)

USING THE THEORY OF HYDRAULIC CIRCUITS IN SIMULATING THERMAL

where T is the transpose symbol, and D, S are the diag onal matrices

D = diag{ d1 , d 2 ,..., d n } , S = diag{ s1, s 2,..., s n}. (2) The matrix A is obtained from the matrix A by deleting any row (usually the last mth row). The row is deleted to make the number of equations equal to the number of unknown variables. Eliminating the second equation in system (1) we obtain

AD = Q; A T P + H = SDD,

pm = pm* , H = H + Am pm* ; Am is the last column of the

matrix AT . The determination of flow distribution in the HC is brought to solving μ equations with μ unknown variables, μ = m + n – 1. As is shown in [1], the solution of system (3) exists and is a unique one. It is usually taken that the pressure drop across each branch with the number υ ∈ [1, 2,..., n] as a function of the flowrate through the branch and its hydraulic resis tance is given by yυ + hυ = sυ dυ dυ ; yυ = pi − p j ,

(4)

where i and j are the input and output nodes of the branch υ. Equations (4) are called the closing rela tions. They form the third vector equation in system (1). In [1], more general closing relations are used: β−1 yυ = sυ ( dυ ) dυ , where the exponent β is the same for any υ. For gases, the following closing relations are sometimes used: yυ + hυ = sυ dυ dυ; yυ = pi2 − p 2j ,

(5)

provided that they are satisfied in all sections of the HC. In this case, Kirchhoff’s first law remains unchanged. Kirchhoff’s second law is replaced by its analog: the sum of squared pressure drops over a closed loop is equal to zero [4]. The general derivations of closing relations for some cases of the flow of com pressible and incompressible fluids can be found in [5]. The method of nodal pressures is one of more attrac tive methods for solving the problem of flow distribu tion in the paths of power installations. The process of searching for the solution by means of successive approximations in the method of nodal pressures con sists in determining the corrections for the pressures at the nodes, which are determined through calculating the discrepancies of flowrates. SPECIFIC FEATURES OF THE PROBLEM CONSIDERED The specific feature of calculating flow distribution in the paths of power installations lies in the closing relations that interconnect the pressures at the HC THERMAL ENGINEERING

Vol. 58

nodes, flowrates through the branches, and the ther mophysical properties of the flow. The expressions describing these dependences, which are various in form and frequently nonlinear in nature, can look as follows [6–7]: pressure drop in a turbine stage 2

No. 13

2011

2

2

(

2

2

)

−2

pi = pi +1 + di +1 pi,0 − pi +1,0 di +1,0Ti /Ti0,

(6)

and steam flowrate through a control valve

(3)

 = ( p1, p2,..., pr ), r = m – 1, because where the vector P the last component of the vector P is usually known:

1095

{

d vlv =

}

0.5 2 2k ⎡− 0.09 p 2 + 1.09 p p − p 2 ⎤ . (7) i i + i i + 1 1 ⎦ 1 − ε cr k + 1 ⎣ In [8], a scheme is proposed for constructing a net work model of flow distribution in the paths of power installations with variable parameters. The authors developed a modified version of the method of nodal pressures by introducing a generalized dependence for the pressure difference as a function of coolant flow rate, and the use of this version makes it possible to take into account the specific features pertinent to the motion of flows in the path elements the pressures in which vary according to nonlinear laws. The algorithm for calculating the corrections in the method of nodal pressures implies the need to deter mine the interconnection between the discrepancies Δ P and ΔD. To do so, the equations of system (3) are linearized, after which the following system of alge braic equations is obtained: A ⎞ ⎛ ΔP ⎞ 0 ⎛ 0 ⎛ ⎞ (8) ⎜ J ( P ) −2SD ⎟ ⎜ ΔD ⎟ = − ⎜ w ( P , D ) ⎟, 0 0⎠⎝ ⎝ ⎝ 2 0 0 ⎠ ⎠ where J ( P0 ) is the Jacobian matrix calculated at the point P0, w2 ( D0, P0 ) = F (P0 ) − SDD ; F(P) is the pressure drop law (4) and (5). Multiplying the second matrix row of system (8) by the multiplier − A [2SD0 ]−1 and sub tracting it from the first row, we obtain the system of the order r (9) C ΔP = b ,

σ vlv

where C = A [2SD0 ]−1 J ( P0 ) ; b = − A [SD0 ]−1 w2 ( D0, P0 ), and we solve it using some technique, e.g., the Gauss method. After that, we calculate Δ D from the formula −1 ΔD = [2SD0 ] ⎡⎣w2 ( D0, P0 ) + J ( P0 ) ΔP ⎤⎦ .

(10)

The use of the described procedure makes it possi ble to take into account closing relations (6) and (7) for (3), which are nonlinear with respect to pressure, while remaining within the framework of the initial method of nodal pressures. Calculations of the circuits of the actual equipment used in power units and steam generators carried out under the conditions of strong disturbances applied to the hydraulic circuit parame ters, such as full closing of a control valve and opera tion of a safety valve, have demonstrated that the com putations carried out using the proposed algorithm for

1096

TAIROV et al.

numbering, these are the last m − r + 1 equations), as a result of which the following system of equations is obtained:

ΔP, MPa 1.35

T AD = Q, A P = Y ; Y + H = S (P, D)U ( P ) DD ,

1.25 65

85

where U ( P ) = diag{u1 ( p) ,..., un ( p)}; uk ( P ) =

QF, %

Effect of heat release QF in the furnace of a PK24 boiler on the pressure loss ΔP in the steamgenerating surfaces at a constant flowrate of coolant.

calculating flow distribution feature fast convergence and good stability. The second specific feature relating to the standard problem is due to the dependence of hydraulic resis tances on the parameters of medium sυ = s( pυ , dυ , xυ ) ,

(11)

where xυ is the steam quality. The approximations of the dependence characteriz ing hydraulic resistance (11) in equipment elements con taining a boiling coolant have the following form [9]: s=

λ ⎛1 + x ψ ⎛ ρ ' − 1⎞⎞, ⎜ ⎜ ⎟⎟ 2 gF ρ ' ⎝ ⎝ ρ '' ⎠⎠

(12)

where ρ ', ρ '' are the densities of water and steam in saturated state, and ψ is the coefficient that takes into account the influence of flow structure on the hydrau lic resistance; this coefficient is a function of steam quality, pressure, and flowrate. The coefficient ψ as a function of the parameters of twophase flow ρ, d, p, and x can be determined from the approximations given in [9]. An analysis of such approximations showed that there are discontinuities in the space of the parameters ρ, d, p, and x. A generalization of what was said above leads us to consideration of an HC for which its own laws of pres sure drop across each circuit with the number υ ∈ [1,2,..., n] are allowed in the general case, in the fol lowing form: 2

fυ ( pυ ) + hυ = s( pυ , dυ , xυ )dυ .

(13) Dependence (13) does not satisfy the requirements of the general HC theory for closing relations [1] because it is not a continuously differentiated one, which adds difficulty to the use of Newton’s method. CALCULATION OF FLOW DISTRIBUTION IN THE CHANNELS OF POWER INSTALLATIONS The problem described above is solved using the double iteration method. We followed the following rule for equating the number of equations to that of unknown variables: the rows with the numbers of the known components of the vector of pressures are dis carded in the subsystem AD = Q (with the selected

S ( P, D ) = diag{ s1( p, d),..., s n( p, d)}.

(14)

pi − p j + hυ ; fυ ( pυ ) + hυ

One of the vector equations contained in equalities (14) can be eliminated, after which we obtain:  T (15) AD = Q; A P + H = S (P, D)DD ,  where S (P, D) = S (P, D)U ( P ) . System (15) consists of r + n equations with r + n unknown variables: n compo nents of the vector D and r components of the vector P, which form the vector P = ( p1, p2,..., pr ). Thus, an ana

log of system (3) is obtained that corresponds to the method of nodal pressures in the HC theory [1, 3], the use of which makes it possible to take into account the specific features of the pressure drop law, including linear and nonlinear ones, all in a single calculation algorithm.  System (16) contains the multiplier S (P, D) , which depends on the closing relations determining the hydraulic pressure in the parts of an HC, also in the form of strongly nonlinear and implicit dependences for sections containing a boiling coolant. Therefore, we will specify in system (15) concrete values of the  vectors S (P, D) and P = P0 in the matrix D = D0 (owing to the stable behavior of the applied method of nodal pressures, the initial approximation P0, D0 may deviate considerably from the sought solution) and organize the iteration process: Pl +1 = Pl + ΔPl и Dl +1 = Dl + ΔDl , l = 0,1,...,

(16)

where the corrections to the flowrates ΔDl and pres sures Δ Pl are determined in solving the system A 0 ⎛ 0 ⎞ ⎛ Δ Pl ⎞ ⎛ ⎞ =−⎜   ⎟⎜ ⎜ T ⎟ ⎟. ⎝ w2(Pl , D l ) ⎠ ⎝ A − 2S (Pl , Dl )Dl ⎠ ⎝ Δ Dl ⎠

(17)

where

 Dl = diag{ d1, d 2,..., d n}, w2(Pl , D l ) =  T A Pl − S (Pl , D l )Dl D. The iteration process continues until the following conditions are satisfied:

Δ Pl ≥ ε p, Δ Dl ≥ ε d ; w2(Pl , D l ) ≥ ε , where ε p, ε d , ε are some specified numbers. The norms of vectors Δ Pl ; Δ Dl ;, w2(Pl , D l ) are understood to mean the maximal moduli of their components. The system of equations (17) takes into account that the last m − r components of the vector y are fixed (specified) and their increments are equal to zero. To simplify the calculations, we multiply the second THERMAL ENGINEERING

Vol. 58

No. 13

2011

USING THE THEORY OF HYDRAULIC CIRCUITS IN SIMULATING THERMAL

1097

Checking the convergence behavior of the modified method of nodal pressures No. of cal Initial pressure culation discrepancy 1 2 3

Pressure discrepancies, MPa, calculated in the iterations according to Newton’s method 1st iteration

2nd iteration

3rd iteration

4th iteration

5th iteration

6th iteration

–1.50818 –1.12487 –1.12382

0.23892 0.12850 0.11501

0.01998 0.00649 0.00604

–0.00149 –0.00114 –0.00173

0.00074 0.00027 0.00026

0.00001 0.00000 0.00000

1.73190 –7.42098 –13.61440

matrix row in system (17) by the multiplier   −1 − A[2S (Pl , Dl )Dl ] and subtract it from the first row. As a result, the system of the order r in the discrepancies Δ Pl is obtained:   −1 T A[2S (Pl , Dl )Dl ] A Δ Pl =   − A[2S (Pl , Dl )Dl ]−1w2(Pl , D l ).

(18)

We will use the solution of this system in calculating ΔDl by the formula  ΔDl = [2S (Pl , D l )]−1[w2 ( Pl , D l ) + A T ΔPl ].

(19) The described scheme is a version of Newton’s method with the difference that system (17) does not include the derivatives of elements of the matrix  S (P, D) . Owing to this circumstance, it becomes much easier to implement numerical algorithms, because there is no need to approximate the derivatives of non linear and implicit (in the case of describing the prop erties of a twophase flow) expressions. With the values  of the matrix S (P, D) calculated in the preceding itera tion, we obtain the classic scheme of calculating flow distribution using the method of nodal pressures. As is shown in [1], such a system has a unique solution, and Newton’s method converges. With the calculation scheme described above, the entire set of branches in the steam–water path of power installations is described in the models by the following single closing relation:

fυ ( z, w ) = c0 z + c1w + c2 z 2 + c3zw + c4w 2 , = sυ (z, w, dυ ,..) dυ dυ

(20)

where w ( z, w) is some known function. By specifying the appropriate coefficients c0–c4 in (20) it is possible to describe the flow of water in tubes, passage of steam through turbine compartments (6), passage of medium through control valves (7), and flow of steam–liquid mixture in heating surfaces or tubes tak ing into account the dependence of hydraulic resis tance [9] on the flow parameters. The calculation scheme was used for determining, in a stepwise manner, flow distribution in the dynamic models of drumtype and oncethrough boilers, tur bine unit, and power unit as a whole. The use of this scheme in constructing mathematical support for a power unit training simulator [10] made it possible to THERMAL ENGINEERING

Vol. 58

No. 13

2011

take into account the change of pressure losses con nected with the value of steam quality in the steam generating heating surfaces of a boiler unit (see the fig ure) at different levels of heat load. For checking the convergence behavior of the cal culation scheme, we carried out test calculations of flow distribution in a boiler unit described by an HC containing 86 branches with specifying the initial approximation P0, D0 essentially different from the sought solution. The results from calculations of pres sure discrepancies at one of the HC nodes are given in the table. It can be seen that convergence is achieved already by the sixth iteration even if the initial approx imation is specified with an essential error, which tes tifies that the proposed method for calculating flow distribution features high stability. CONCLUSIONS A further development of the method of nodal pressures is presented, which is used as the main tech nique for solving flow distribution problems in the hydraulic circuit theory taking into account the spe cific features of various nonlinear and in a number of cases also implicit closing relations describing the motion of coolant in the steam–water paths of power installations. The computational algorithm developed on the basis of the proposed procedure has shown sta ble performance and high speed of operation when applied to example calculations of the process circuits of a steam boiler and the entire power unit. The described procedure for calculating flow distribution can be implemented as part of full highspeed mathe matical models of power installations used at thermal power stations. ACKNOWLEDGMENTS This work was supported by the Russian Founda tion for Basic Research, grant no. 090800201. REFERENCES 1. A. P. Merenkov and V. Ya. Khasilev, Theory of Hydraulic Circuits (Nauka, Moscow, 1985) [in Russian]. 2. O. K. Shashkov and V. O. Shashkov, “Calculating Vari able Operating Conditions of Thermal Power Stations with SteamTurbine Installations Using the Method of

1098

3. 4. 5. 6. 7.

TAIROV et al.

Thermohydraulic Circuits,” Therm. Eng. 51, 325 (2004). E. V. Sennova and V. G. Sidler, Mathematical Simula tion and Optimization of Heat Supply Systems (Nauka, Novosibirsk, 1987) [in Russian]. K. R. AidaZade, “Computational Problems in Hydraulic Networks,” Zh. Vychisl. Mat., Matem. Fiz. 29 (2), 184–193 (1989). O. A. Balyshev and E. A. Tairov, Analysis of Transients and Steady Processes in Pipeline Systems (Nauka, Novosibirsk, 1998) [in Russian]. F. A. Vul’man and N. S. Khor’kov, Thermal Calcula tions of Thermal Power Installations on a Computer (Energiya, Moscow, 1971) [in Russian]. V. A. Ivanov, Operating Modes of Large Steam Turbine Units (Energiya, Moscow, 1971) [in Russian].

8. E. A. Tairov, V. F. Chistyakov, and I. V. Karaulova, “Application of a Network Model for Calculating Flow Distribution in the Paths of Power Installations,” Izv. Ross. Akad. Nauk, Energetika, No. 3, 105–113 (2003). 9. D. F. Peterson, L. V. Saminskaya, and V. B. Khabenskii, “On the Development of a Method for Carrying Out Hydraulic Calculation of a SubcriticalPressure Once Through Boiler on a Computer,” Trudy TsKTI, Issue 98, 60–70 (1969). 10. A. A. Levin and E. A. Tairov, “Applying the Method of Nodal Pressures for Calculating Flow Distribution in the Paths of Power Installations,” in Proceedings of the Seminar of Higher Schools in Siberia and Far East “Problems of Thermal Physics and Thermal Power Engi neering,” Irkutsk, 2008, pp. 46–49.

THERMAL ENGINEERING

Vol. 58

No. 13

2011