Estimating friction errors in MOC analyses of unsteady pipe flows

Computers & Fluids 36 (2007) 1235–1246 www.elsevier.com/locate/compfluid

Estimating friction errors in MOC analyses of unsteady pipe flows Masashi Shimada a, Jim Brown b

b,*

, Alan Vardy

b

a Graduate School of Life Sciences, University of Tsukuba, Japan Division of Civil Engineering, University of Dundee, Nethergate, Dundee DD1 4HN, UK

Received 15 March 2005; received in revised form 1 July 2006; accepted 12 November 2006 Available online 9 February 2007

Abstract Errors arising from the numerical treatment of friction in unsteady flows in small pipe networks are assessed for fixed-grid, method of characteristics analyses (MOC) with no interpolation. Although the results of the study are targeted at general unsteady flows, the underpinning analytical development is based on the behaviour of standing waves. This enables quantitative conclusions to be drawn about the sensitivity of MOC solutions to frequency. The development is undertaken first for a single pipe and then for networks, with and without a loop. The performance of MOC solutions is modelled with the aid of polynomial transfer matrices that are analogous to transfer function matrices available directly from partial differential equations of water hammer. It is found that, in a single pipe, the numerical errors tend to increase damping at all frequencies, but that, in networks, either increased or decreased damping may occur. The errors place both higher and lower limits on the frequencies that can propagate along a pipe in a numerical analysis. This contrasts with the physical case where only a lower limit exists. The work presented is part of a wider project with the long-term aim of automating the selection of numerical grid sizes in MOC analyses of unsteady flows. Ó 2007 Elsevier Ltd. All rights reserved.

1. Introduction The analysis of unsteady flows in pipe systems is a challenging subject that involves at least three key sources of potential error. First, it is necessary to choose simplified geometrical and physical characteristics to make possible a manageable theoretical treatment. Second, there is the representation of the simplified system by analytical expressions. Third there is the solution of the equations, usually by numerical means. The first of these sources commonly includes considerable simplification of boundary conditions, especially pumps and valves, etc., as well as neglecting physical phenomena such as cavitation. The second can include neglecting two- and three-dimensional effects and making simplifying assumptions about physical phenomena – assuming uniform wavespeed and neglecting the dependence of skin friction on rates of change of the *

Corresponding author. Tel.: +44 0 1382 384505; fax: +44 0 1382 384816. E-mail address: [email protected] (J. Brown). 0045-7930/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2006.11.005

flow state, for instance. The third includes the method of representing analytical equations by numerical counterparts and then the method of integrating the numerical equations. This paper addresses one aspect of the latter topic, namely errors arising from the integration of skin friction terms in analyses based on the method of characteristics (MOC). The work is part of a wider study with the (admittedly ambitious) long-term target of making it possible to automate the selection of numerical grids. When using software packages for analyses of unsteady flows in pipe networks, analysts usually need to pay attention to both fluid dynamics and potential sources of numerical error. In an ideal world, these tasks would be separated so that engineers could focus on the fluid dynamics and the software would provide error limitation automatically. This paper is a small step towards that aim. It targets rational predictions of errors caused by the numerical integration of friction terms in MOC analyses of flows in networks. Shimada and Vardy [16] made a preliminary approach to this issue, but were forced to resort to approximate methods to allow for junctions and other boundaries.

1236

M. Shimada et al. / Computers & Fluids 36 (2007) 1235–1246

Nomenclature A B c D EA FD Fp fD Gp g H h i j Kk k L M MOC m N n PA Q q R s t x Dt X Dx

cross-sectional area (m2)  c=gA wave speed (m/s) eigen-decomposition of FD network exact augmented transfer matrix polynomial transfer matrix over a grid length polynomial transfer matrix for a pipe with reference index p Darcy friction factor exact transfer matrix for a pipe gravitational acceleration (m/s2) pressure head (m) perturbed pressure head (m) complex number, i2 ¼ 1 grid length index angle corresponding to one grid length of undamped kth mode mode number pipe length (m) modulus of the complex amplitude w method of characteristics number of time steps in the period of a numerical oscillation number of grid lengths in a pipe number of time increments (t ¼ nDt) network polynomial augmented transfer matrix volumetric flow rate (m3/s) perturbed volumetric flow rate (m3/s) resistance factor defined in Eq. (2) complex frequency (s1) time coordinate (s) distance along a pipe (m) time interval in numerical grid (s) matrix of eigenvectors distance interval in numerical grid (m)

Herein, attention is focussed on the fixed grid method of characteristics [25] and, for clarity and simplicity, there is no interpolation. In principle, errors due to interpolation can be addressed in a similar manner, but it is helpful to understand the individual contributions of friction and interpolation (and other causes of error) before combining them in a general analysis. In a similar way, the authors have addressed the influence of interpolation in the absence of friction [17,18]. Numerical errors due to interpolation have been investigated by Wiggert and Sundquist [22] and by Goldberg and Wylie [6] for space-line and time-line interpolation, respectively. In both cases, Fourier methods were used to deduce amplitude and phase errors in inviscid flows in single pipes. These studies provide a solid basis for the development of grid assessment methodologies in more complex flows – e.g. including friction and pipe junctions. Fourier (or sim-

Y Zc

vector of network physical variables characteristic impedance

Greeks b; b1 ; b2 bL c c2 v d s h K k1;2 g r m U w x X

coefficients – see Eqs. (5) and (6)  Nb propagation constant defined in Eq. (27) phase ratio amplitude performance ratio of MOC eigenvector component (Appendix A) numerical phase angle for one time step wave length (m) eigenvalues of the polynomial transfer matrix FD real part of complex frequency (Eq. (28)) wave number (m1) kinematic viscosity (m2/s) exponent for one time step – see Eq. (30) complex amplifying factor angular frequency (s1) exact phase angle for one time step

Subscripts bdy boundary j grid point index (j ¼ 1 : N þ 1) k mode number lam laminar flow n nth time increment ðt ¼ nDtÞ ‘; r pipe indices turb turbulent flow u, d up-, down-stream Superscript steady flow *

ilar) analyses were also used by Maudsley [13], Lai [10] and Yang and Hsu [27]. Another approach is to assess numerical error during post-processing (i.e. assess the errors after performing the simulation). This is a powerful method of gaining general understanding. For example, Chaudhry and Hussaini [3] showed that second-order accurate schemes offer relatively little value with Courant numbers close to unity, but are of increasing benefit with decreasing Courant number. Sibetheros et al. [19] used spline interpolation in a frictionless flow and found significant improvement over linear interpolation and over second-order finite difference schemes. In the absence of friction, the Courant number is a sufficient criterion to describe stability. In the presence of friction, however, a second condition must usually be satisfied. Explanations of the need for the second condition and

M. Shimada et al. / Computers & Fluids 36 (2007) 1235–1246

practical methods of allowing for it are given by Holloway and Chaudhry [8], Anderson et al. [1] and Wylie [26]. Karney [9] proposed a method of assessing and adjusting grid size errors during unsteady flow analyses, thereby avoiding unnecessary computations during periods where changes occur only slowly. He achieved this by monitoring rates of change of internal energy and kinetic energy. Dominance of the former implies a need for water-hammer methods of analysis whereas dominance of the latter implies that rigid-column methods are adequate. A similar outcome was achieved by Murray [14], although less generally. Ghidaoui and Karney [4] and Ghidaoui et al. [5] demonstrated the underlying effects of using time-line and space-line interpolation. They examined the use of wave speed adjustment techniques as an alternative to interpolation. Selek et al. [15] have compared numerical computations with measurements of water hammer pressure transients in a penstock. They used fixed grid MOC with no interpolation and with space-line interpolation, and variable grid MOC; friction was incorporated explicitly. The computed results showed good agreement with the measurements, the performance being best with variable grid and poorest with simple MOC. None of the above papers enables a priori prediction of numerical errors in unsteady flows with friction in pipe networks. The purpose of this paper is to assess errors in predicted amplitude and phase exactly for transient MOC analyses in the presence of friction, but with no interpolation. That is, the difference between the actual numerical solution and the corresponding analytical solution of the same equations is obtained exactly. In the following development, attention is directed firstly at flows in single pipes and then at flows in networks. The chosen methodology uses a linearised approach based on standing wave theory and so it is able to provide predictions only for a finite number of discrete frequencies. Nevertheless, these frequencies cover the whole range of reasonable interest for users of one-dimensional methods. The MOC predictions are obtained using polynomial transfer matrices and the errors are determined by comparison with exact solutions of the governing analytical equations. In this way, accurate estimates are made of the influence of discretization on the accuracy of the numerical method. Amplitude and phase errors are presented for simple networks with and without a loop. Throughout the paper, attention is focused on the behaviour of wave propagation along pipes and the implications of this for grid size selection. In principle, the cause of the unsteadiness is irrelevant to this purpose even though it is of fundamental importance in any particular flow application. Accordingly, we model only the simple case of passive boundary conditions such as a prescribed constant pressure or a prescribed constant flow rate. The prescribed unsteady state conditions already exist at time zero and we study the influence of the numerical grid structure on the subsequent decay. It is expected, but not pro-

1237

ven, that the results (or, at least, the trends) will be valid for any passive boundary conditions. They may also be valid in the presence of active boundary conditions, but this is not certain. Furthermore, special attention would need to be paid to grid size selection even in the absence of any form of damping (numerical or physical) because the interactions at boundaries will depend on discretization of the waveforms at the boundaries. 2. Stability and error analysis of MOC equations in a single pipe Fig. 1 illustrates paths of MOC lines in adjacent reaches of a regular ðx  tÞ grid. Values of the flow parameters are to be obtained at the point A based on previously calculated values at the points L and R. The gradients of the characteristic lines C þ and C  are equal to the wave speed c. To avoid interpolation, the grid dimensions Dx and Dt are chosen to satisfy c ¼ Dx=Dt and so the Courant number is equal to unity. A uniform grid is considered so that Dx ¼ L=N where L is the length of the pipe and N is the number of grid lengths. In general network analysis, it is not usually possible to choose a Courant number of unity in every pipe, but this problem is avoided herein by the judicious choice of pipe lengths and wave speeds. The need for interpolation in general networks is considered in a previous paper [18]. The following analysis assesses the response of a steady flow with head H  and discharge Q , to perturbations h, q in the head and discharge. Thus the total head and discharge are H ¼ H  þ h and Q ¼ Q þ q. For this case, the MOC equations can be written in terms of the perturbed flow parameters as (see [25]) Z A Z A Z A Cþ : dh þ B dq ¼  Rq dx ð1aÞ L

L

L

and C : 

Z

A

dh þ B

Z

R

A

dq ¼ 

Z

A

ð1bÞ

Rq dx

R

R

in which B ¼ c=gA, g is the gravitational acceleration, A is the pipe cross-sectional area and R is a resistance coefficient per unit length of pipe. The particular forms of R for laminar and turbulent flow (suitably linearized) are taken as

A

n

C+

Δt

C-

L

j -1

R Δx

j

Δx

Fig. 1. Definition sketch for characteristics lines.

n -1 j +1

1238

Rlam ¼

M. Shimada et al. / Computers & Fluids 36 (2007) 1235–1246

32m gAD2

and

Rturb ¼

fD jQ j ; gDA2

ð2a; bÞ

respectively. In both cases, this implies that the skin friction is approximated by quasi-steady expressions even though much is now known about the dependence of frictional resistance on historical flow conditions as well as upon the instantaneous velocity [20,21,28]. For the purposes of this paper, it is considered more useful to focus on quasisteady approximations because they remain so widely used. The need to linearise the friction term also merits comment. In particular, in the case of turbulent flow, it restricts the validity of the analysis to small perturbations (although some behaviour predicted for larger perturbations may be instructive qualitatively). No such restriction is implied in the laminar flow case (Eq. (2a)) provided, of course, that the amplitudes are not so large as to cause breakdown of the laminar structure of the flow. The left-hand sides of Eqs. (1a) and (1b) are exact differentials, but the right-hand sides require integration. Since the variation of q is not known a priori, numerical integration is necessary and the outcome depends upon the particular approximation chosen. Many alternatives are possible (see for example [24]) and each has its advantages and disadvantages. For the purposes of this paper, a trapezoidal rule approximation is adopted. That is, the mean perturbed skin friction along a characteristic line in any time step is regarded as a simple average of its values at the beginning and end of the time step. This approximation, which further limits both laminar and turbulent analyses to small perturbations, leads to the finite difference equations: C þ : hj;n  hj1;n1 þ Bðqj;n  qj1;n1 Þ þ ðR=2Þðqj;n þ qj1;n1 ÞDx ¼ 0

ð3aÞ

and C  : hj;n þ hjþ1;n1 þ Bðqj;n  qjþ1;n1 Þ þ ðR=2Þðqj;n þ qjþ1;n1 ÞDx ¼ 0

ð3bÞ

in which the indices ðj; nÞ identify, respectively, the spatial and temporal positions in the fixed rectangular grid (see Fig. 1). It is convenient to re-express these equations as C þ : hj;n  hj1;n1 þ Bðb1 qj;n  b2 qj1;n1 Þ ¼ 0

ð4aÞ

and C  : hj;n þ hjþ1;n1 þ Bðb1 qj;n  b2 qjþ1;n1 Þ ¼ 0;

ð4bÞ

where b1 ¼ 1 þ b;

b2 ¼ 1  b;

ð5a; bÞ

and ðbÞlam ¼

16mDx cD2

and

ðbÞturb ¼

fD jQ jDx : 2cDA

tice, the second of these requirements is more demanding than the first. The parameter ðbÞturb is equivalent to the parameter w used by Wylie [26] who accepted it as an adequate representation of the frictional effect and found that it should be less than unity for stability. In contrast, values of less than 0.1 (perhaps much less) are generally needed to constrain numerical damping adequately, but there is no formal method of quantifying the required value. One purpose of the work begun in this paper is to reduce the need for empiricism in the choice of b. The turbulent resistance parameter ðbÞturb can be interpreted in several ways as may be seen in the following rearrangement: ðbÞturb ¼

fD Dx jQ j 12 fD ðDx=DÞðQ =AÞjQ =Aj=g ¼ cðQ =AÞ=g 2D cA

ð7Þ

The first identity can be interrelated as the product of a dimensionless grid length friction parameter and a Mach number. The second identity can be interpreted as the quotient of the grid length friction head and the Joukowski head for the quasi-steady flow Q . Corresponding interpretations can be demonstrated for the laminar case. As seen in Eq. (6), b is proportional to the grid length Dx. For future reference, it is useful to have an equivalent parameter that is independent of the numerical grid. Accordingly, the parameter bL ¼ N b is defined using the pipe length in place of the grid length in Eqs. (6a,b). In the turbulent case, this has the form ðbL Þturb ¼ N ðbÞturb ¼

fD L jQ j ; 2D cA

ð8Þ

which is equivalent to the parameter R studied by Liou [11,12] and by Wylie [26] in the context of fast transients. Wylie [26] considered that this parameter should be supplemented by one depending on the rate of change of, say, velocity. In total, there are N pairs of MOC equations of the form of Eqs. (4a,b) for each pipe. Together with boundary conditions at the ends of the pipe, there are 2ðN þ 1Þ equations to be solved simultaneously to give the perturbed head and flowrate at each of the N þ 1 grid points. With linear friction, these equations can, in principle, be solved for any linear boundary condition. For the turbulent case, however, exact solutions are possible only when the mean flow component is steady. This condition is satisfied for perturbed flows with the form of standing waves and, in practice, it is necessary to choose cases with either a node or an anti-node at each end of the pipe, i.e. hbdy ¼ 0 or qbdy ¼ 0. As an example, the boundary conditions for a reservoir-pipe-valve system could take the form h1;n ¼ 0

and

qN þ1;n ¼ 0:

ð9Þ

ð6a; bÞ

The parameter b, which is always greater than zero, has an important role in the behaviour of the numerical equations. Analysts must ensure that it is sufficiently small to ensure stability and to constrain numerical damping and, in prac-

2.1. Stability analysis Consider a liquid-filled pipe, unbounded in its axial direction and subjected to a disturbance with wave number

M. Shimada et al. / Computers & Fluids 36 (2007) 1235–1246

r ¼ 2p=K, where K is the wave length. The response of the MOC equations to the disturbance has the form hj;n ¼ h0 wn expðirjDxÞ and qj;n ¼ q0 wn expðirjDxÞ;

ð10a; bÞ

where h0 and q0 are the initial magnitudes of the perturbation head and flow rate. The parameter w is, in general, complex and describes the temporal evolution of the original wave, i.e. its changes in amplitude and phase. Substituting Eq. (10) into Eq. (4) leads to a stability polynomial, namely b1 w2  2w cosðrDxÞ þ b2 ¼ 0:

ð11Þ

This equation is valid for any wave number r and the following identities may usefully be noted 2p 2p Dx ¼ cDt ¼ xDt; ð12Þ K K where X is the angular change corresponding to one grid length or time step and x is the angular frequency. In the case of standing wave perturbations in a finite pipe, however, only certain discrete frequencies are possible. To distinguish these from the general case, an integer subscript k is applied to wave properties to denote the kth standing wave mode. The frequencies of the damped physical oscillation are then defined by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Xk  ðK 2k  b2 Þ; ð13Þ

X  rDx ¼

where K k is given, for odd and even standing wave modes, respectively, by   1 p or Kk ¼ k  2 N kp Kk ¼ for k ¼ 1; 2; 3 . . . ; N ; ð14a; bÞ N and is the angular advance corresponding to one grid length of the undamped kth mode. Eq. (13) implies that, when K k < b, no oscillatory solution is possible. This corresponds to the region to the left of the straight line in Fig. 2. 1.5 over-damped (physical & numerical)

1.25

critical damping (physical)

damped oscillations (physical) over-damped (numerical)

1 critical damping (numerical)

β

0.75 mode 4 N=32,16,8,4

damped oscillations (physical & numerical)

The solutions of Eq. (11) fall into one of three categories, each corresponding to a particular type of numerical behaviour, namely damped, critically damped and (socalled) over-damped. These cases are discussed independently in the following paragraphs. Case 1: 0 < b < sin Xk Case 1 is the only category in which standing waves persist in the numerical solution. It is obtained for sufficiently small values of b and for a limited range of frequencies. The standing wave patterns continue indefinitely, but their amplitudes reduce monotonically. For this case, Eq. (11) has complex conjugate roots, namely qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos Xk þ i sin2 Xk  b2 w¼  Meihk and 1þb qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos Xk  i sin2 Xk  b2  Meið2phk Þ ; w¼ 1þb ð15a; bÞ with modulus sffiffiffiffiffiffiffiffiffiffiffi 1b M  jwj ¼ 1þb

ð16Þ

and argument 0qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 sin2 Xk  b2 1 @ A: hk  arg w ¼ tan cos Xk

ð17Þ

The modulus and argument of w determine the amplitude and frequency of the oscillation described by the numerical solution. Since b < 1 (in Case 1), it follows from Eq. (16) that the predicted oscillations will decay, which is consistent with physical reality for standing waves in a viscous medium in the absence of an energy source. Case 2: 0 < b ¼ sin Xk The second category of solution is the special case of critical numerical damping, namely just sufficient to prevent the oscillating behaviour implied by standing waves. That is, the perturbed conditions die away without oscillation (although one sign change is possible with certain initial conditions). For this case, Eq. (11) has two equal roots, namely   sffiffiffiffiffiffiffiffiffiffiffi cos Xk   ¼ 1  b  M: jwj ¼  ð18Þ 1þb 1þb

A

0.5

modes 1-8 N=8

0.25

0

1239

0

0.5

1

1.5

2

2.5

3

Kk

Fig. 2. Regions of oscillating and non-oscillating flow.

3.5

Clearly M < 1 and the amplitude decays as M n , where n denotes the number of time steps since the chosen initial condition. The critical damping condition is displayed as the curved line in Fig. 2. It is the locus of the condition b ¼ sin Xk . The region below the curve corresponds to the damped oscillations discussed in Case 1.

1240

M. Shimada et al. / Computers & Fluids 36 (2007) 1235–1246

Case 3: 0 < sin Xk < b The third category of possible solutions occupies the region in Fig. 2 outside the critical numerical damping curve. Much of this region corresponds to large values of b where the skin friction terms are so great that competent analysts would know instinctively that the MOC analysis is unsafe. Even with small values of b, however, there are two regions where the numerical solution will not exhibit continuing oscillation even though exact analytical solutions would do so.This case has two sub-categories. The simpler of these represents an over-damped system. Once again, the perturbed conditions die away with, at most, only one change of sign. In this case, however, they die away more quickly than in the case with critical damping. The other sub-category is a non-physical oscillation with a time period equal to two time steps. For both sub-categories, the roots of the complex amplitude are real, namely qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cos Xk  b2  sin2 Xk w¼ ð19Þ 1þb with jwj < 1 for b > 0. Depending on the value of cos(Xk ), both roots may be positive, both may be negative or they may be of opposite sign. Positive roots correspond to the first sub-category and the amplitude decreases with M n where n is the number of time steps since the prescribed initial condition. Negative roots correspond to the second sub-category, namely a non-physical, decaying oscillation whose frequency is a function of the numerical grid, not of the flow conditions.

"

h1;n

#

 ¼

q1;n

F 11

F 12

F 21

F 22

"

hN þ1;n qN þ1;n

#

"  Fp

hN þ1;n

#

qN þ1;n

in which the outermost grid points are denoted by j ¼ 1 and j ¼ N þ 1, respectively. For any particular physical case, the elements of Fp could be calculated numerically by evaluating individual terms in the matrix FD and taking successive products to N obtain ðFD Þ . Alternatively, they can be found analytically using an eigen-decomposition of FD . It is shown in Appendix A that the pipe transfer matrix Fp depends on the eigenvalues of FD and on the ratio of its off-diagonal terms. The resulting expressions are 1 F 11 ¼ F 22 ¼ ðkN1 þ kN 1 Þ; 2 s F 12 ¼ ðkN1  kN 1 Þ 2

ð23aÞ ð23bÞ

and F 21 ¼

1 N ðk  kN 1 Þ; 2s 1

ð23cÞ

in which s and the eigenvalue k1 are defined in Appendix A. Since k1 depends on w2 , it can be seen that the elements of the pipe transfer matrix are polynomials of degree 2N. The roots of these polynomials determine the numerical response to a given set of end conditions. That is, they describe the amplitude decay rate and the frequency error. 2.3. Exact standing wave solutions An exact analysis of free oscillations – e.g. [25] – leads to the counterpart of Eq. (22) in the form 

hð0; tÞ qð0; tÞ



 ¼

cosh cL

Z c sinh cL

ðsinh cLÞ=Z c

cosh cL



hðL; tÞ qðL; tÞ

qj;n ¼ wqj;n1 :

ð20a; bÞ

By combining these expressions with Eq. (4), the system variables at adjacent spatial grid positions j and (j þ 1) can be shown to satisfy "

hj;n qj;n

#



1 ¼ 2w

"

b1 w2 þ b2 Bðb21 w2  b22 Þ ðw2 1Þ B

b1 w2 þ b2

 Gp

hðL; tÞ qðL; tÞ

 ;

ð24Þ

The preceding section deals with temporal responses at any particular location in a pipe. Attention now turns to relationships between conditions at adjacent spatial points at any particular time. First, Eqs. (10a,b) are used to relate conditions at successive times at a typical location j, namely and





2.2. Polynomial transfer matrix

hj;n ¼ whj;n1

ð22Þ

;

#"

hjþ1;n qjþ1;n

#

"  FD

hjþ1;n qjþ1;n

# ;

ð21Þ

where the square matrix FD is a transfer matrix for the grid length. For a finite pipe with N equal grid lengths, there are N such equations and they can be combined to give a transfer matrix Fp ¼ ðFD ÞN for the whole pipe, namely

where c and Zc are, respectively, the propagation constant and the characteristic impedance, defined by c¼

s pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 2bL =ðsL=cÞ c

and

Zc ¼

ðc=sÞ ; ðgA=c2 Þ

ð25a; bÞ

in which s is the complex frequency. The pipe transfer function matrix Gp is the exact solution counterpart of the polynomial transfer matrix Fp. As an example, for oscillatory solutions in a reservoirpipe-valve system with hð0; tÞ ¼ 0 and qðL; tÞ ¼ 0 (i.e. odd mode shapes), the propagation constant is purely imaginary (i.e. c ¼ ic2 ) and is determined from coshðcLÞ ¼ coshðic2 LÞ ¼ 0:

ð26Þ

The propagation constant for this case, with k ¼ 1; 2; 3; . . ., and referring to Eqs. (12) and (14), is     1 p 2p Kk KkN ¼ ic2 ¼ i ; ð27Þ ¼i c¼i k ¼ irk ¼ i 2 L Kk L Dx

M. Shimada et al. / Computers & Fluids 36 (2007) 1235–1246

and the complex frequency is v" ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi#ffi u  2  2  2 u bL 1 p bL t þi k  s¼ 2 L=c ðL=cÞ ðL=cÞ  g þ ixk :

ð28Þ

For the kth mode, the angle Xk and the amplitude exponent U for a time interval equivalent to one time step in the numerical solution can be obtained by multiplying Eq. (28) by Dt  L=Nc and using Eq. (27) and bL ¼ N b to yield vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi u" 2   2 # u 2 1 p bL Xk ¼ xk Dt ¼ t k   2 N N qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ ðK 2k  b2 Þ ð29Þ and U ¼ gDt ¼ bL =N ¼ b:

ð30Þ

From Eq. (29), it is seen that oscillatory solutions exist only for modes k for which b < K k . The critical physical damping condition is given by b ¼ K k and the amplitude decreases in time without oscillations when b > K k . By inspection, these categories are closely analogous to those described in Section 2.1 for the numerical analysis, but the limits are different and, of course, there is no sub-category corresponding to non-physical oscillations. In Fig. 2, the physical critical damping condition is seen as the straight line b ¼ K k . To its right, the conditions reflect damped oscillations. To its left, there are no oscillations. 3. Comparison of MOC prediction with exact analysis for a single pipe MOC predictions for standing waves in a single pipe are now compared with exact solutions of the perturbation (linearised) equations. The purpose is to gain information about the error sensitivity as a function of frequency both for amplitude and phase. In nearly all unsteady flows of practical interest, a wide spectrum of frequencies will be present, whereas, in this analysis, only a discrete number of frequencies can be investigated explicitly – namely those for which standing waves exist with a node or anti-node at each end of the pipe. This is a somewhat artificial constraint, however, because exactly the same analysis would apply if the finite length of pipe under consideration were part of an infinite pipe. In the latter case, any physically reasonable wavelength can be considered and the trends are monotonic (see, for instance [6,22]). The intermediate frequencies could alternatively be studied by using the fixed length pipe in series with one or more other pipes in the manner described in the following section on pipe networks. Eq. (15) and (16) show that the amplitude of the MOC oscillatory solution changes by a factor of M in one time step. Similarly, Eq. (30) shows that the exact oscillation

1241

changes by a factor of expðgDtÞ ¼ expðbÞ ¼ expðbL =N Þ. The ratio of these two factors, designated herein as a performance ratio of the MOC solution, is sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffi 1  bL =N 1b 1 þ b 1þb =N M L ¼ : ð31Þ d¼ U¼ eb ebL =N e In this expression, bL is a property of the pipe (and of the underlying steady flow in the case of turbulent flow) and b is additionally dependent upon the numerical grid (see (6) and (8)). The ratio d tends towards unity when b tends to zero, which implies successive reductions of the viscous resistance and/or of the grid size. In general, therefore, the accuracy of the solution is controlled by the value of b, not by the grid size or the system properties alone. For example, performance ratios of d ¼ 0:99 and d ¼ 0:95 correspond to b ¼ 0:3053 and b ¼ 0:5059 respectively. As discussed in the paragraph following Eqs. (6), b is much smaller than unity in practical engineering applications. Accordingly, the behaviour of Eq. (31) may be approximated by series expansion as d  1  b3 =3 . . . ;

ð32Þ

which shows that the numerical amplitude-error grows very slowly when b is small. These equations also show that the ratio d is less than unity, and so the numerical waves decay more rapidly than the true waves. From Eq. (31), it is seen that the amplitude performance ratio d is independent of the mode number k. In Fig. 2, a constant amplitude error ratio implies a horizontal line spanning the region beneath the curve denoting critical numerical damping. The independence of the performance ratio d on the mode number k arises because d is defined per time step in Eq. (31). The performance per wave period is poorer at high frequencies than at lower frequencies. The ratio of the numerical and analytical frequencies is obtained from Eqs. (17) and (29) as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sin2 Xk  b2 1 tan hk cos X qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k v¼ : ð33Þ ¼ Xk 2 K  b2 k

This ratio is found to be sensitive to mode number, grid size and friction parameter. The fractional error v  1 is usually small, but can be either positive (i.e. numerically increased frequency, decreased period) or negative (numerically decreased frequency, increased period). Numerical experiments show that (v  1) is negative for sufficiently low modes and sufficiently coarse grids. It tends to become less negative, and can become positive, with increasing mode number and/or with finer grids (i.e. larger N). 3.1. Consequences for numerical analysis It is normal for numerical solutions to differ quantitatively from analytical solutions. For flows of the type

1242

M. Shimada et al. / Computers & Fluids 36 (2007) 1235–1246

considered herein, however, Fig. 2 shows that strong qualitative differences may also exist. That is, there are large regions where the numerical solution exhibits no oscillations even though these exist in practice. Consider, for instance, the line of open circles corresponding to modes 1–8 for a particular value of b and N ¼ 8. Mode 1 is over-damped both physically and numerically. Mode 8 is over-damped numerically although physically it should be a damped oscillation. Modes 2–7 are oscillatory, but experience both physical and numerical damping. The influence of the grid size is illustrated by the filled triangles in Fig. 2. These correspond to mode 4 with N ¼ 4, 8, 16, and 32. The points lie on a straight line through the origin. For the particular mode shown, the numerical solution is over-damped when N ¼ 4, whereas it exhibits damped oscillations with larger N. The improved accuracy per time step with increasing N more than compensates for the implied increase in the number of time steps required for each complete wave period. For completeness, attention is drawn to the point denoted by ‘‘A’’ and depicted by a filled square in Fig. 2. This point corresponds to mode 2 simulated with N ¼ 8 and an appropriate value of b. The analytical solution corresponds to damped oscillations whereas the numerical solution implies over-damped behaviour. The only difference between this point and the open circle at the same value of K k in Fig. 2 is a reduction of 12.5% in the pipe diameter. This has the effect of increasing the friction parameter b by about 40%, which, in this instance, is sufficient to change the predicted numerical behaviour. It is instructive to note that the over-damping cannot be corrected by decreasing the grid size because, as discussed in the preceding paragraph, the locus of the filled square for the same mode would be a straight line through the origin.

its grid point N p þ 1) is connected to the left hand ends of pipes l and r (i.e. their grid points 1). At a typical instant t ¼ nDt, continuity of mass and (for simplicity) equality of heads at the junction can be expressed as

4. Pipe networks

respectively. Eq. (35) may be written as

Although the preceding discussion was focussed on a single pipe, it has strong implications for the analysis of flows in networks. The results are not directly transferable, however, because the conditions at internal junctions are not, in general, nodes or anti-nodes. To explore behaviour in a network, it is necessary to develop transfer matrices for the complete system. This is done by assembling a network transfer matrix from the transfer matrices for individual pipes, supplemented by appropriate equations describing flows and heads at internal junctions and external boundaries. Book and Watson [2] describe alternative schemes for defining transfer matrices. Herein, the chosen methodology is a simple extension of the method described above for single pipes. The process is the same as that used by the authors [18] in a related study of interpolation with different Courant numbers in the various pipes. A network transfer matrix can be assembled from transfer matrices for the individual pipes together with appropriate matching conditions at the internal junctions. Suppose that the right hand end of a typical pipe p (i.e.

PA Y ¼ 0;

ðqp ÞN p þ1;n ¼ ðq‘ Þ1;n þ ðqr Þ1;n

ð34aÞ

and ðhp ÞN p þ1;n ¼ ðh‘ Þ1;n ¼ ðhr Þ1;n :

ð34bÞ

In the special case of two pipes p and l connected in series, the terms involving hr and qr are omitted from the above equations. Likewise, when more than three pipes meet at a junction, additional terms are required. At an external boundary, one side of one of Eqs. (34a) and (34b) will be a prescribed quantity. The complete set of equations comprising transfer matrices, internal junction conditions and external boundary conditions may be expressed schematically in matrix form as 32 Y 3 2 u1 I F1 0 0 0  6Y 7 7 6 0 7 6 d1 0 I F2 0    76 6 76 Y 7 6 7 6 0 u2 7 0 0   7 76 6 6 ð35Þ 76 Yd2 7 6 Jn Cnds 7¼0 76 6 7 7 6 . .. 7 .. 76 6 . 7 56 4 . . 4 . 5 .. Bdy Cnds . in which ‘‘Jn Cnds’’ are the junction conditions, ‘‘Bdy Cnds’’ are the boundary conditions and the physical variables at, say, the upstream end of pipe 1 and the downstream end of pipe 2 are represented as     hu1 hd2 Yu1 ¼ and Yd2 ¼ ; ð36a; bÞ qu1 qd2 ð37Þ

where PA (a function of w) is the network polynomial transfer matrix, the subscript A indicating that it is augmented by the junction and boundary conditions. The vector Y comprises a list of pairs of head and discharge values at the ends of the pipes. A non-trivial solution of Eq. (37) is found using the determinant condition jPA j ¼ 0;

ð38Þ

which has the form of a polynomial in w. The roots of this polynomial equation define the modes of response of the network. A similar procedure using the free oscillation transfer function matrix Gp (Eq. (24)) leads to equivalent augmented equations EA Y ¼ 0

ð39Þ

where EA is the exact network transfer matrix augmented by the junction and external boundary conditions. The solution of

M. Shimada et al. / Computers & Fluids 36 (2007) 1235–1246

jEA j ¼ 0

ð40Þ

yields values of the complex frequency s, which has real and imaginary parts giving damping and frequency properties for the modes of oscillation. 5. Numerical examples The network equations are now applied to two specific examples using nodes and antinodes at external boundaries. The topologies of the two cases are illustrated in Fig. 3. One case is a tree-like branched system with 5 pipes and the other is a looped system with 6 pipes. To obtain suitable initial conditions, each system is first solved for the underlying steady state flows, thereby enabling the values of Q to be determined. Then, polynomial and exact transfer matrices for the complete system are obtained in the manner described in Section 4. In principle, many alternative formulations are possible, depending upon the choices of nodes and anti-nodes at the various external boundaries. In the particular examples considered herein, constant head conditions are specified at locations 1 and 2 and constant flow conditions are assumed at location 6. The results for the branched and branched/loop networks are discussed in Sections 5.1 and 5.2, respectively. For completeness, amplitude and phase errors have also been deduced from actual MOC simulations. These are not strictly necessary, but they confirm the validity of the polynomial predictions and they illustrate a difficulty faced by all analysts when attempting to interpret numerical predictions. For the MOC analysis it is necessary to provide valid initial conditions of head and flow rate at the grid positions in each pipe for each mode of oscillation of the network. The appropriate complex amplifying factor w is known from the solution of Eq. (38). Using this in the grid length polynomial transfer matrix Eq. (21) and stepping progressively inwards along the pipes from the boundaries enables the initial values at the interior grid points to be generated. The MOC analysis generates values of head and discharge at discrete intervals that do not, in general, include the exact instants at which maxima, minima and null conditions occur. To infer these instants, it is necessary to make estimates by interpolation between values obtained in successive time steps. Herein, a simple linear interpola-

1

pipe 3

pipe 6

5

4 2

tion has been used. Greater accuracy might be achievable using other forms of interpolation for this post-processing activity, but this is considered unnecessary for present purposes. Frequencies can be assessed from instants at which the (interpolated) head or discharge crosses zero-head or zero-discharge state. Amplitudes are deduced at instants when the discharge-time graph crosses the zero-discharge state. 5.1. Branched network with no loop Data used for the branched system example are given in Table 1. In the initial steady flow, the heads at locations 1, 2 and 6 are 23.5793 m, 23.008 m and 0 m, respectively, giving a discharge of 1.0472 m3/s at location 6. In the perturbed flow, the head has nodes at locations 1 and 2 and an anti-node at location 6. This implies that the perturbed discharge has anti-nodes at locations 1 and 2 and a node at location 6. Fig. 4a shows the amplitude error (ratio of error to true amplitude) per wave period for three values of the time step. In addition, the error per time step is shown for the longest and shortest of the time steps. In each case, all possible standing wave modes are shown. The lowest mode is determined by the pipe length, but the highest mode increases with decreasing grid length. The errors are quite small, even with the coarsest grid – for which two pipes have only one grid length, two others have only two grid lengths and the longest pipe has three grid lengths. In contrast with the single pipe case, the errors are positive for some modes (i.e. increased decay rate), but negative for others (decreased decay rate). The errors increase in absolute value as the mode increases, and decrease with decreasing time step. The errors per wave at Dt ¼ 0:8 s are noticeably larger than per time step but the difference is relatively smaller for higher modes as there are fewer time steps in one wave. At Dt ¼ 0:2 s the difference are hardly noticeable on the scale of the figure.

Table 1 Pipe and flow parameters for network examples Case

Pipe

Length (m)

Diameter (m)

Darcy f

Wavespeed (m/s)

Steady flow (m3/s)

Branch

1 2 3 4 5

2400 1600 800 1600 800

1 1 0.7 0.8 0.6

0.03 0.03 0.03 0.03 0.02

1000 1000 1000 1000 1000

0.5269 0.5203 0.5269 0.5203 1.0472

Loop

1 2 3 4 5 6

2400 1600 800 1600 800 3200

1 1 0.7 0.8 0.6 0.9

0.03 0.03 0.03 0.03 0.02 0.03

1000 1000 1000 1000 1000 1000

0.5254 0.4818 0.5277 0.5195 1.0472 0.0377

pipe 1 3

pipe 5

6

pipe 4

pipe 2

Fig. 3. Topology of branched and looped network systems.

1243

The networks used in these examples have been chosen solely for the purposes of illustrating the numerical method. They are not representative of any particular application.

1244

M. Shimada et al. / Computers & Fluids 36 (2007) 1235–1246 0.012

0.008

m

2.0E-4

Dt=0.8s 1wave Dt=0.4s 1wave Dt=0.2s 1 wave Dt=0.2s 1step Dt=0.8s 1step

0.004

Dt=0.8s Dt=0.4s Dt=0.2s

1.6E-4

1.2E-4

0

8.0E-5

-0.004

4.0E-5

-0.008

0.0E+0 0

5

10 k,

s

15

20

-1

0

5

10 k,

s

15

20

-1

Fig. 4. Influence of grid size on numerical accuracy – branched system. (a) Amplitude errors per wave period ðm ¼ 2p=hk Þ, and time step ðm ¼ 1Þ, (b) frequency errors.

When the simulations for any particular mode are repeated with different time steps (not shown in Fig. 4), it is found that the predicted decay rate increases strongly with Dt. This is because b itself increases with the grid size. For the lowest oscillation modes, comparisons of individual numerical data for particular modes show that the 2 amplitude error grows approximately as ðDtÞ . For higher modes, the error depends even more strongly on the time step. This behaviour is highlighted in Fig. 5, which shows the same data expressed as absolute errors and presented as a function of the number of time steps per period for each mode. If the error dependence on Dt were linear, this figure would show no dependence on the grid size because any changes in the error per time step would be exactly counterbalanced by corresponding changes in the number of time steps per wave period. Fig. 6a gives a comparison between the numerical decay per wave period predicted by the polynomial method and the decay obtained in actual MOC simulations. As indicated above, the two values should be identical, but practical difficulties confuse the interpretation of actual simulations, especially at higher frequencies. The close correlation between the two methods of determining the numerical errors confirms the validity of the polynomial methodology, thereby justifying its use in grid-selection processes. The particular results shown in Fig. 6a were 0.012

Dt=0.8s Dt=0.4s Dt=0.2s

|

m

|

0.008

0.004

0 1

10

100

T k / Δt

Fig. 5. Dependence of absolute amplitude errors per wave period on the number of time steps in the wave period ðm ¼ 2p=hk Þ – branched system.

obtained using the smallest grid size considered in Fig. 4a (Dt ¼ 0:2 s). Similar agreement is obtained with the coarser grids. Fig. 4b shows the errors in numerically predicted frequencies for each possible mode using the various grid sizes. The errors are very small, even with the coarsest grid and all correspond to a contraction of the period of the oscillation – i.e. to an increase in frequency. The trend lines in the figure are included for visual impact only. They show that the errors, albeit very small, increase approximately 2 with ðDtÞ , which is consistent with the behaviour found for single pipes. Fig. 6b compares predictions of the numerical phase errors. The open squares show values deduced by post-processing of actual numerical simulations and the filled circles shows values obtained using the polynomial approach. The agreement is generally satisfactory except at two points (x  4 and x  9) at which the post-processing methodology is seen to be weak. 5.2. Branched network with a loop Table 1 shows the data used for the 6-pipe system with a simple loop. There are some differences from the data used for the previous example. The initial steady flow heads at locations 1 and 2 are both equal to 23.8385 m and the initial head at location 6 is 0 m. These imply a discharge of 1.0472 m3/s at location 6 that remains constant throughout the subsequent unsteady flow. In common with the branched system example, the perturbed head has nodes at locations 1 and 2 and an anti-node at location 6 whereas the perturbed discharge has anti-nodes at locations 1 and 2 and a node at location 6. Inferred errors in the numerical predictions of the amplitude decay per wave period are shown in Fig. 7a. In common with the results for the simple branched system, the errors are small and may be either positive or negative. That is, although damping is positive in all cases, the absolute value is sometimes greater and sometimes smaller than the analytical value. In common with the branched system, for any particular mode, the errors vary approximately as ðDtÞ2 . Also, for any particular mode of oscillation in the

M. Shimada et al. / Computers & Fluids 36 (2007) 1235–1246

a

0.16

b 1.5E-05

Polynomial Simulation 0.12

1.0E-05

0.08

5.0E-06

0.04

0.0E+00

0

1245

Polynomial Simulation

-5.0E-06 0

2

4

6 k,

8

s

10

12

0

2

4

-1

6 k,

8

10

12

s -1

Fig. 6. Comparison of polynomial predictions and actual simulations ðDt ¼ 0:2 s) – branched system. (a) Amplitude decay in one wave period, (b) frequency change.

a

b

0.02

Dt=0.8s Dt=0.4s Dt=0.2s

0.015

1.6E-4

Dt=0.8s Dt=0.4s Dt=0.2s

1.2E-4

m

0.01 0.005

8.0E-5 0 -0.005

4.0E-5

-0.01

0.0E+0

-0.015 0

5

10 k,

s

15

20

-1

0

2

4 k,

6

s

8

-1

Fig. 7. Influence of grid size on numerical accuracy – looped system. (a) Amplitude errors per wave period ðm ¼ 2p=hk Þ, (b) frequency errors.

network, the amplitude decay rate increases with the time step – because it increases with b, which itself increases with Dx and hence with Dt. The fractional frequency error of the polynomial predictions is shown for three time steps in Fig. 7b. For the smallest time step, the errors are very small and increase little over the range of frequencies shown. Once again, the general tendency is for the error in phase of a given mode of oscillation to increase approximately in proportion to the square of the time step. Although not shown, the level of agreement between the polynomial predictions and actual simulations is comparable to that shown in Fig. 6 for the branched system with no loop. 6. Conclusions 1. The decay of unsteadiness in flows along pipes with friction has been studied through the superposition of small-amplitude perturbed flows on underlying steady flows. The perturbed flows have the form of standing waves with nodes or anti-nodes at all external boundaries, but not at internal boundaries (junctions). 2. It has been shown that the locus of the critical damping curve predicted by MOC analyses for a single pipe can differ markedly from its analytical counterpart. One effect is to extend the true region of monotonic damping at low frequencies. Another effect is to introduce a non-

physical region of monotonic damping at high frequencies. Oscillatory damping is numerically possible only between these regions and it always exceeds the true physical damping. 3. Polynomial transfer matrices have been introduced to study the corresponding behaviour in pipe networks. The network solutions provide amplitude and frequency information that characterise the performance of full MOC analyses. Comparisons of these data with corresponding data from exact analytical solutions enable the performance of alternative numerical grids to be studied. 4. In contrast with the behaviour found for single pipes, amplitude errors in networks can be either positive or negative. For typical numerical predictions, even with quite coarse grids, amplitude errors are quite small and phase errors are very small. Both types of error increase approximately with the square of the grid size. 5. Predictions from the polynomial approach have been compared with equivalent predictions obtained by post-processing results of direct simulations using an MOC analysis. The two sets of predictions agree closely for a range of grid sizes and for a range of oscillation frequencies. Isolated cases of poor agreement are attributable to difficulties in post-processing that face any analyst seeking to infer amplitudes and frequencies from discrete numerical data.

1246

M. Shimada et al. / Computers & Fluids 36 (2007) 1235–1246

Acknowledgments The authors express gratitude for financial assistance from the Japanese Society for the Promotion of Science (Grant-in-Aid No. 14360140) and from the European Commission’s Fifth Framework Growth Programme (ref: Surge-Net, G1RT-CT-2002-05069). The authors are solely responsible for the content of the paper, which does not necessarily reflect the opinion of JSPS or the EC. Appendix A. Nth power of the polynomial transfer matrix The following method of expressing the Nth power of a polynomial matrix summarises fuller presentations in, say, Golub [7] and Wilkinson [23]. Let the elements of the polynomial transfer matrix FD be denoted by f11 , f12 , f21 , f22 . The eigenvalues are pffiffiffiffiffiffiffiffiffiffiffi ðA:1Þ k1;2 ¼ f11  f12 f21 and it is easily verified that k1 k2 ¼ 1: ðA:2Þ The matrix of eigenvectors is   s s X¼ ðA:3Þ 1 1 where pffiffiffiffiffiffiffiffiffiffiffiffiffiffi s ¼ f12 =f21 : The eigen-decomposition of matrix FD is thus   k1 0 1 D ¼ X FD X ¼ 0 k2 and hence its Nth power is " # 0 kN1 N 1 N : D ¼ X FD X ¼ 0 kN2 Finally, XDN X1

ðA:4Þ

ðA:5Þ

ðA:6Þ

" # N N N N k þ k sðk  k Þ 1 1 2 1 2 ¼ FND ¼  Fp : 2 ðkN1  kN2 Þ=s kN1 þ kN2 ðA:7Þ

References [1] Anderson A, Arfaie M, Sandoval-Pena R, Suwan K. Pipe friction in waterhammer computation. In: Proc XXIV IAHR congress, Madrid, Spain, vol 2; 1991. p. 23–30. [2] Book WJ, Watson C. Alternatives in the generation of time domain models of fluid lines using frequency domain techniques. Math Comput Simul 2000;53:353–65. [3] Chaudhry MH, Hussaini MY. Second-order accurate explicit finitedifference schemes for waterhammer analysis. J Fluids Eng ASME 1985;107:523–9. [4] Ghidaoui MS, Karney B. Equivalent differential equations in fixedgrid characteristics method. J Hydraulic Eng ASCE 1994;120: 1159–75.

[5] Ghidaoui MS, Karney B, McInnis DA. Energy estimates for discretization errors in water hammer problems. J Hydraulic Eng ASCE 1998;124:384–93. [6] Goldberg DE, Wylie EB. Characteristics method using time-line interpolations. J Hydraulic Eng ASCE 1983;109:670–83. [7] Golub GH. Matrix computations. 2nd ed. Baltimore (MD): The John Hopkins University Press; 1990. p. 369. [8] Holloway MB, Chaudhry MH. Stability and accuracy of waterhammer analysis. Adv Water Resour 1985;8:121–8. [9] Karney B. Energy relations in transient closed-conduit flow. J Hydraulic Eng ASCE 1990;116:1180–96. [10] Lai C. Comprehensive method of characteristics models for flow simulation. J Hydraulic Eng ASCE 1989;114:1074–97. [11] Liou CP. Maximum pressure head due to linear valve closure. J Fluids Eng Trans ASME 1991;113:643–7. [12] Liou CP. Adequacy of waterhammer equations for leak detection in petroleum products pipelines. In: Proc international conference on unsteady flow and fluid transients, Durham, UK, 29 September–1 October, HR Wallingford; 1992. p. 355–69. [13] Maudsley D. Errors in the simulation of pressure transients in a hydraulic system. Trans Inst Meas Contr 1984;6:7–12. [14] Murray SJ. Transient multi-regime pipe flow – a case study. In: Proc international conference on unsteady flow and fluid transients, Durham, UK, 29 September–1 October, HR Wallingford; 1992. p. 121–36. [15] Selek B, Kirkgo¨z MS, Selek Z. Comparison of computed water hammer pressures with test results for the C ¸ atalan power plant in Turkey. Canad J Civil Eng 2004;31:78–85. [16] Shimada M, Vardy AE. Towards automatic error limitation in MOC analyses. Proceedings of the 8th international conference on pressure surges, The Hague, Netherlands, 12–14 April. Cranfield, UK: BHR Group; 2000. p. 183–204. [17] Shimada M, Vardy AE, Brown JMB. Numerical accuracy – manual vs automatic assessment. Proceedings of the 9th international conference on pressure surges, Chester, UK, 24–26 March. Cranfield, UK: BHR Group; 2004. p. 511–26. [18] Shimada M, Brown JMB, Leslie D, Vardy AE. Time-line interpolation errors in pipe networks. J Hydraulic Eng ASCE 2006;132:294–306. [19] Sibetheros IA, Holley ER, Branski JM. Spline interpolations for waterhammer analysis. J Hydraulic Eng ASCE 1991;117:1332–51. [20] Vardy AE, Brown JMB. Transient turbulent friction in smooth-pipe flows. J Sound Vibr 2003;259:1011–36. [21] Vardy AE, Brown JMB. Transient turbulent friction in fully rough pipe flows. J Sound Vibr 2004;270:233–57. [22] Wiggert DC, Sundquist MJ. Fixed-grid characteristics for pipeline transients. J Hydraulic Division ASCE 1977;103:1403–16. [23] Wilkinson JH. The algebraic eigenvalue problem. Oxford: Oxford University Press; 1965. p. 33. [24] Wylie EB. Advances in the use of MOC in unsteady pipeline flow. Proceedings of the 4th international conference on pressure surges, Bath, UK, 21–23 September. Cranfield, UK: BHR Group; 1983. p. 27–37. [25] Wylie EB, Streeter VL, Suo L. Fluid transients in systems. Englewood Cliffs (NJ): Prentice Hall; 1993. [26] Wylie EB. Unsteady internal flows - dimensionless numbers and time constants. In: Proceedings of the 7th international conference on pressure surges and fluid transients in pipelines and open channels, Harrogate, UK, 16–18 April. Cranfield, UK: BHR Group, Mechanical Engineering Publications; 1996. p. 283–8. [27] Yang JC, Hsu EL. Time-line interpolation for solution of the dispersion equation. J. Hydraulic Res IAHR 1990;28:503–20. [28] Zielke W. Frequency-dependent friction in transient pipe flow. J Basic Eng 1968;90:109–15.