- Email: [email protected]
Quasi-steady friction in transient polytropic flow
Pergsmon
Computers& Figi& Vol. 26, No. 8, pp. 793-809, 1997 0 1997 ElsevierScience Ltd. All rights reserved Printed in Great Britain 00457930/97 $17.00 t 0.00
PII: s0045-7930(97)00029-7
QUASI-STEADY
FRICTION ALAN
IN TRANSIENT
VARDYt
and ZHIQIU
POLYTROPIC
FLOW
PAN
Civil Engineering Department, University of Dundee, Dundee DDl 4HN, U.K. (Received 24 October 1996; in revised farm 24 October 1997) Abstract-An accurate method of integrating quasi-steady friction in transient polytropic flows is introduced and is justified both physically and analytically. When applied to subsonic flows, it is found to be especially informative at higher Mach numbers where existing methods of integration can lead to serious errors in mass conservation, The method is used herein to assess the accuracy of widely used finite difference representations of quasi-steady friction, attention focusing prin~pally on solutions obtained using the method of characte~sti~. The new method is exact in the special case of steady flows. Of the approximate methods considered, the popular approximation V,/V, _ AJAXis found to perform well. A modification introduced by Karney and McInnis (Journal of Hydology Engineering, ACSE, 1992, 118(7), 1014-1030 [I]) improves the correlation for d’ownstream characteristic lines, but has the opposite effect for upstream characteristics. A simple cubic. scheme, loosely related to a more complex one used successfully by Murray (Proceedings of the International Conference on &steady Flow and Fhdd Transients, 29 September-l October. Balkema, Amsterdam, 1992, pp. 143-157 [2]), illustrates difficulties in higher-order schemes. The new method is not exact in transient Rows, but it gives good accuracy in slowly varying flows and it is plausible in rapidly varying flows. Using it as a benchmark, the approximation k’,lV, _ 4,iAr is again shown to perform well. All of the schemes considered are unreliable when the flow is about to become cho:ked within the region of integration. 0 1997 Elsevier Science Ltd.
NOMENCLATURE A C
C+ c s LHS M MOC n P RHS t V x it Ax Y P f Subscripts A av B dn L, R I r-At up W
cross-sectional area [m*] polytropic wavespeed [=&p/p)] [m s-‘1 characteristic lines in the direction of flow characteristic lines against the direction of flow pipe diameter [m] skin friction efficient [ET,,,/(O.Sp yi Vi)] cross-sectional perimeter [m] left-hand side of MOC equation polytropic Mach number [=V/c] method of characteristics polytropic index pressure [Pa] right-hand side of MOC equation time coordinate [s] average veiocity in cross-section [m P’] distance coordinate [m] elevation [m] time interval (e.g. tA-tL) [s] distance interval (e.g. xA-xL) [m] ratio of principal specific heat capacities mass density kg me31 shear stress [Pa] solution point average value intermediate point on MOC line downstream end of sofution region base points on MOC lines solution time previous calculation time upstream end of solution region wall
~Corresponding author. 793
A. Vardy and Z. Pan
194
1. INTRODUCTION
The method of characteristics (MOC) is the most widely used numerical scheme for simulating transient liquid flows in pipes, especially when there are numerous boundaries. It is capable of yielding accurate predictions for a wide range of flow types provided that limitations on grid sizes are recognized [3]. The picture is less satisfactory for flows in which there are significant variations in wavespeed or where particle velocities are important in comparison with wavespeeds. Common examples include gas flows in long pipelines and free-surface flows in channels. In such cases, the usual methods of expressing MOC equations numerically can lead to serious errors in mass conservation [2,4]. This paper addresses one cause of inaccuracies in MOC when Mach numbers or Froude numbers are significant, namely the integration of terms such as skin friction. One purpose is to introduce a new, highly accurate method of integrating these terms. A more important purpose is to use the new method to assess the validity of simpler, commonly used methods. 1.1. Sources of error in MOC In their simplest form, MOC equations can be integrated exactly. This is so when terms on the right-hand sides of equations such as Equations (4) and (5) herein are zero. It can be possible to obtain exact solutions for simple boundary conditions. In practical applications, the equations must be integrated numerically because of complications from terms describing friction, area change, heat transfer, etc., and from difficulties in locating characteristic lines-especially in fixed-grid methods of solution. Murray [2] set out the three main causes of inaccuracy, namely (i) Location of characteristic lines-the (analytical) properties of the “left-hand side” (LHS) terms in the MOC equations apply only in precisely defined directions in the (x-t) plane. In practice, it is impossible to deduce these directions exactly at all necessary locations. (ii) Interpolation-usually, parameter values must be deduced at locations other than those at which they are notionally known. The resulting interpolation processes involve numerical approximation. (iii) Integration of RHS terms-the numerical integration of the right-hand side (RHS) terms involves approximations when the variations of parameters (e.g. velocity) cannot be expressed analytically. All three of these difficulties can exist in a general method of analysis and all can cause significant errors. Their relative importance is difficult to assess a priori and it is sometimes difficult to distinguish between them. For the purposes of this paper, however, attention is focused exclusively on the third item. Errors due to the first two causes are minimized (often eliminated) by careful selection of the test cases. Where possible, analytical solutions are used to provide exact parameter values at the limits of the region of integration. For simplicity of presentation, the discussion is restricted to polytropic gas flows in horizontal pipes of constant area. In principle, however, the method can be used in many other applications and also in other numerical schemes, e.g. finite differences and finite elements. 2. METHOD
OF CHARACTERISTICS
(MOC)
The equations of continuity and momentum for one-dimensional velocity V along a duct of non-uniform area A may be written as
flow of a perfect gas at a
(1) and
(2) in which p is the mass density, p is the absolute pressure, r, is the wall shear stress, &,, is the length of the duct perimeter, g is the gravitational acceleration and dz/dx is the elevation gradi-
195
Transient polytropic flow
Fig. 1. Characteristic lines
in the (x-r) plane(C’, dx/dt =
V + c; C,
dx/dr = V- c).
ent. Both dz/dx and dA/dx are expressed as an ordinary derivatives to indicate that changes in time are excluded. herein. For the purposes of this paper, all flows are assumed to be polytropic. The effective wave speed c is assumed to satisfy c2 = Y P where n is the polytropic index. The following development is applicable when the flow is not isothermal, that is n # 1. The special case of isothermal flow is presented in Appendix A. Equations (1) and (2) are partial differential equations in two independent variables (x and t). In MOC, they are recast as ordinary differential equations, valid only when specific relationships exist between x and t, namely: dx when -= dt
2
V+c:
dc
dV
(4)
and dx when - = V-c: dt
2 de dV dz --dt=+g;i;;+--7dx. n-ldt
zW&
cVdA (5)
PA
In MOC, the co:nditions dx/dt = I/k c are used to construct characteristic lines in (x-t) space. In Fig. 1, for example, the line LA is chosen so that its gradient at each point satisfies dx/ dt = V + c. As a consequence, Equation (4) is valid at all points on the line and so, in principle, it may be integrated from L to A. Similarly, Equation (5) may be integrated from R to A. 2.1. Integration qf LHS terms The integration of the left-hand side (LHS) terms in Equations (4) and (5) with respect to time t is trivial and the results are analytically exact. Along LA for example:
--$g+g
dt= 1
-& (CA - CL)+ (VA -
VL).
(6)
No assumptions are required in this integration. In particular, no information is required about the manner in w’hich c and V vary along LA (except that they do so continuously). This is the source of the popularity of MOC. It is not always possible to express the equations in a form where the integration can be performed exactly, but it is often possible to do so. In other cases, it m.ight be necessary to use numerical approximations in, for example, coefficients of the term dc/dr.
796
A. Vardy and Z. Pan
0
0.2
0.4
0.6
0.8
1
k/L Fig. 2. Steady polytropic flow in a constant-area duct cfL/O = 50, n = 1.4, pdn = 100 kPa, Tdn= 300K).
2.2. Integration of RHS terms
The integration of the RHS terms in Equations (4) and (5) is intrinsically more difficult than that of the LHS terms because of the existence of expressions such as cV and VIVI. To integrate these exactly would require a knowledge of their variations in time, which is not possible a priori. In practice, therefore, the integration must be undertaken numerically. There are two circumstances in which the need for numerical evaluation does not have a serious adverse influence on the final solution. The first is when the magnitudes of the RHS terms are much smaller than the magnitudes of the LHS terms. This is quite common in highly transient, low Mach number flows. It is less common in higher Mach number flows, partly because transient events disperse more readily and partly because velocity-dependent terms tend to be much larger. The second circumstance is when the true variation of, say, VIVJ is well represented by the numerical approximation. This is quite common in regions of low Mach number flow remote from transient disturbances. In such regions, there is little variation in velocity on any particular characteristic line and so all reasonable approximations to VIVI yield similar values. At higher Mach numbers, however, considerable variations in velocity can exist, even in a steady flow, so different methods of approximating VIVI can yield significantly different values. One purpose of this paper is to extend the range of cases when the numerical integration can be performed accurately-by removing the restriction to low Mach number flows. This is achieved by providing a good approximation to the true variation of velocity along characteristic lines. 2.3. Qualitative d@erence between LHS and RHS terms The distinction between the LHS and RHS terms in Equations (4) and (5) is emphasized herein by referring to the latter as environmental terms. These describe background conditions that exist even in the absence of waves (i.e. in a steady flow). Mathematically, they may be functions of the dependent parameters, but they must not depend on derivatives of these parameters. Physically, they are more closely related to steady-flow considerations than to unsteady ones. Small waves tend to have little influence on them. In contrast, the LHS terms are strongly dependent on unsteadiness; they can be large even when the RHS terms are zero. 3. STEADY
POLYTROPIC
FLOW
The detailed development begins with the special case of steady flow. This is important in its own right because many unsteady flows ultimately become steady. It is also the condition in
Transient polytropic flow
797
which lack of m,ass conservation is most easily detected. The steady flow case has even greater significance herem, however, because it underpins the methodology adopted for unsteady flows in Section 4.
3.1. Integration between prescribed Mach numbers Figure 2 shows the axial velocity dist~bution in a steady, polytropic flow of air along a horizontal duct of constant area and constant friction coefficient f (~~~~~0.5~~1~1~). Near the upstream end (low Mach numbers), the velocity increases gradually as the reduction in pressure due to friction causes a decrease in density. Near the downstream end (high Mach numbers), the rates of change are much greater because the skin friction stresses are greater. At any position along the duct, the local rate of change of the polytropic Mach number M with distance x satisfies dM S(n+l) dx=D
N3 (l-M*)’
(7)
Using this expression, the integral s,” I/ * dt may be developed along a characteristic line satisfying dx/dt = Y i- c as
JA~x
=
,I_ l+M
dx
R 1-M
D
=f(n
J LCM
(8)
d”.
To integrate the last of these expressions exactly, it is necessary to eliminate the wavespeed c. The method used to do this depends upon the value of the polytropic index n. The particular case of isothermal flow (n = 1) is considered in Appendix A. The procedure for other values of n is as follows. For polytrop.ic flows, the product c@ - ‘)‘w + ‘) is constant. This enables the wavespeed c inside the integral in Equation (7) to be eliminated, giving n-l A
J
ff*&
L
n+l DcAM~
A
1)
I_L
f(n+
=
I-M
2n Mn+l
d&I
Similarly, along the characteristic line RA n-l
J
A
R
V*dt=-
DC&~ f@+l)
A
1+&f
J R --3TdM _
Mn+l
A. Vardy and Z. Pan
798
=-
DQM;+~ f
3.2. Integration between prescribed locations
The preceding analysis enables s;’ V 2 dt and f,” V2 dt to be evaluated between flow sections at which the Mach numbers are known. In practical applications, however, it is more usual to integrate between prescribed locations than between prescribed Mach numbers. In such cases, the Mach number is known at one end of the region, but it is not known a priori at the other end. Fortunately, the determination of the other Mach number is straightforward. Its value is obtained by integrating Equation (7), giving
The existence of the logarithm in this expression complicates the numerical determination of the unknown Mach number because it necessitates an iterative method of solution. This is undesirable for present purposes because the integration is performed many times in the proposed method of analysis for transient flows (see Section 4). A simple way of overcoming this complication is outlined in Appendix B. 3.3. Comparison with finite d@erence approximations Equations (7)-(11) are exact. Inaccuracies can arise when they are used in MOC analysesbecause of inaccurate values of iwL or xA-xL, say (see Section l.l)-but the integrations themselves are exact. They may therefore be used to assess the accuracy of commonly used approximations. The most common ways of integrating environmental terms utilize finite difference approximations such as (along LA): A JL
V] VI dtx( Vi V&“At.
02)
in which At = (tA--tL) and, by implication, the integration is carried out along a characteristic line. Many ways of approximating (VI VI),, are possible. Examples include (V, _ AtlVr _ J), (V,l Vg/) and (V,,l V,,i) where Vay= 1/2(V, _ dt+ k;). The choices are essentially arbitrary. It cannot be known a priori which will give the best result. Two approximations used extensively herein are (VlVl),,AtX
k’,Iv,-AttAt
(13)
and (VI VI)&=
f V;-Ar
+ @“-t -
vi-At))
f v,-AdAt.
(14)
The approximation VJ V, _ ~~1was presented by Wylie [5] and by Vardy and Chan [6] at the same conference. It is simple, easy to code and it does not necessitate an iterative method of solution. Most importantly, it appears to be unconditionally stable [5,7,8]. Equation (14) was introduced by Karney and McInnis [l]. It is a weighted sum of Equation (13) and the simple approximation V, _ &I Vt _ &it. When c = 1.0, it becomes identical to Equation (13). Karney and McInnis tentatively recommended the use of t; = 0.85, This value is implied herein whenever Equation (14) is used. Murray [2] recognized that the integration could be performed more accurately if values of the velocity could be estimated at intermediate positions on the characteristic lines. The inte-
Transient polytropic flow
0.15
I
(v)
799
I
c+
0.1
0.1
3 il
0 a 0.05
0.05
1
,t=L 0
0.1
0
0.2
0.3
0.4
0.5
f_ Fig. 3. Influence of fAx/D on the average value of Y(vl [n = 1.4, Mdn= 0.9, Tdn= 300K; Sl, v, _ & v, _ A,[; S2, v,l v,l; S3, v,l v, _ ~~1;S4, Karney and McInnis; SS, cubic; S6, polytropic (exact)].
gration could then be carried out in, say, four intervals of 1/4At instead of one interval of At. To obtain the intermediate velocities, he first located the x-positions of the 1/4At points, using dx/dt = Vf c at the ends L and A to obtain a cubic x(t). He then deduced the corresponding depths (JJ) in his free-surface flow, using estimated values of dy/dt at the ends L and A to obtain a cubic y(x). Finally, the corresponding velocities were obtained, based on an assumed linear distribution of fIowrate. Murray’s method is not applied herein. In part, this is because his description intentionally left some decisions open (he wanted to propose a methodology rather than define a universally applicable imple:mentation). In addition, it is because the method involves adjustments that give rise to discontinuous behaviour when presented in the manner used herein. Instead of applying Murray’s method directly, a cubic polynomial V(t) is created using (a) V and t at the ends of the characteristic lines and (b) assumed velocity gradients iilL’/& at these locations. The assumed gradients satisfy
A. Vardy and Z. Pan
800
Table 1. Base values for comparisons of CV*dr 0.25 1.4 0.9 300K
Skin friction fAx/D Polytropic index n Downstream Mach number Mdn Downstream temperature Tdn
at=[ax 1 s dt av
av
--
dx
(15)
where [aV/ax], denotes the gradient that would exist in a steady flow with the same velocity and dx/dt is the gradient of the appropriate characteristic line. Murray began with these gradients, but did not use [aV/ax], directly. Instead, he first made a plausible adjustment and then limited the allowable values to avoid (a) inflexions and (b) the implied existence, within the range t - At to t, of velocities outwith the range V! _ At to V,. These adjustments are not made herein, thus simplifying the presentation, but at the expense of failing to constrain undesirable behaviour when the Mach number approaches unity, implying [aV/ax],+oo. Having obtained a cubic polynomial satisfying the velocities and the gradients [aV/ax],, the evaluation of the integral JV’ dt between the end states at, say, L and A is straightforward. In the remainder of this paper, this is referred to simply as the ‘cubic’ scheme. 3.3.1. Influence offsx/d. Figure 3 illustrates the influence of the parameter fAx/D on the accuracy of five methods of approximation (defined in the figure caption). The top two boxes show the average values of VI VI implied by the various schemes along the characteristic lines LA and RA, respectively. The required values of the integrals $ V 2 dt and s,” V 2 dt are the products of (VI VI),, and the time intervals At z (t,-tl _ A, ) shown in the bottom two boxes. The interval tA - tR consistently exceeds tA-tL because the simulated direction of flow is always from L to R. Upstream-moving waves travel more slowly than downstream-moving waves. The middle pair of boxes shows percentage errors associated with three of these schemes. In all figures herein, the polytropic case is shown as a continuous line. In the lower boxes, it coincides with error = 0%. Except for fAx/D, all parameters in Fig. 3 have the base-case values listed in Table 1 (e.g. n = 1.4). Values for the polytropic (exact) case are obtained for the C+ characteristic line LA as follows: (i) given fAx/D and the values of all parameters at the (downstream) end A, evaluate the interval At = tA-tL by integrating dt = dx/( V i- c); (ii) evaluate all parameters at the upstream end L using Equation (11) and polytropic state relationships; (iii) evaluate jf V 2 dt using Equation (9); (iv) divide by At to give (VI Vl)aY. A similar procedure is used for the C characteristic line RA, for which the point A is chosen at the upstream end of the same duct. All results shown in Fig. 3 have been obtained using exact values of all parameters at the various limits of integration. That is, all end states satisfy analytical relationships and the paths of the characteristic lines are exact. This ensures a fair comparison of the various methods of integration, with no complications arising from location or interpolation errors. Schemes Sl and S2 use values of II VI determined at only one end of the duct. They are seen to be significantly less accurate than the other combinations and so they are not considered hereafter. Because of their relevance to Scheme S5 [Equation (14)], however, note that each overestimates along one characteristic, but underestimates along the other. The curves obtained using Equations (13) and (14) (with E = 0.85) are the most accurate of the approximations considered in the figure, but the errors are significant when fAx/D is sufficiently large. Expressed as a percentage of (VI VI),, , however, the error resulting from the use of Equation (13) is smaller along the C characteristic than along the C+ characteristic. The absolute errors in the integral JV” dt are approximately equal along the two characteristics, how-
Transient polytropic flow
0
0.2
0.4
0.6
0.8
801
0
1
0.2
0
0.2
0.4
0.6
0.8
0.4
0.6
0.8
1
M
M
1
-10I 0
0.2
M
0.4
0.6
0.8
1
M
0.15 8
0
0.2
0.4
0.6 M
0.8
1
0
c-
0.2
I
0.4
0.6
0.8
1
M
Fig. 4. Influence of Mach number on the average value of v V’j[fAx/D = 0.25, n = 1.4, Td,,= 300K; S3, V,l V, _ *,I; S4, Karney and McInnis; S5, cubic; S6, polytropic (exact)].
ever, because of the influence of the time interval At. In contrast, Equation (14) gives approximately equal percentage errors along the two characteristics, but the absolute error jV* dt is much larger along the C characteristic than along the C+ characteristic. Thus the improvement reported by Karney and McInnis for transient, low Mach number flows is not obtained in steady, higher Mach number flows. For all steady flows considered herein, it would be useful to choose 6 < 1 along downstream characteristics, but t: > 1 along upstream characteristics. The detailed quantification of this process would merit a full investigation in its own right because the application of this elegant method is so simple. The cubic scheme performs poorly on the downstream characteristic, but very well on the upstream characteristic. These different behaviours are attributable to the different influences of the gradients of the characteristic lines used in Equation (1 S), but the authors have been unable to find a satisfactory way of exploiting this fact. Adjustments such as those proposed by
802
A. Vardy and 2. Pan
1
1.2
1.4
1.6
1.8
2
1.4
1.2
1
30 1
I
1.6
1.8
R
II
30 ,
2o._.. “‘...I”...... ..^ .. ..YL . l+.._._ _........ I
I
c-
I
il______~ ____.~
;_!
-20
(iv)
2
-- -- --,__,__,______\__-_ *- -- ,______.__.__,__.__
3.
1.2
1
1.4
1.6
1.8
2
II
1
1.2
1.4
1.6
1.8
II
0.15 /
1
2
0.15 (vi) c+-
0.1 -. $ i5
0.M -.
01 1
:
: 1.2
:
: 1.4
: n
: 1.6
: 1.8
:
1 2
04 1
:
1.2
:
:
:
1.4
:
:
1.6
;
1.8
2
n
Fig. 5. Influence of the polytropic index on the average value of fl vi @&.x/B= 0.25, k&=0.9, T& = 300K; 53, YtlV, _ &I; 54, Karney and McInnis; S.5, cubic; 56, polytropic (exact)].
Murray [2] have been found to reduce the errors on the downstream characteristic, but their application renders the process cumbersome and the integral does not necessarily increase monotonically with ~&X/D. Moreover, in the authors’ implementations, the same adjustments had a detrimental influence on the integration along the upstream characteristic. 3.3.2. Injhence of Much number. Figure 4 shows the influence of the Mach number when all other parameters have the base-case values listed in Table 1. The abscissa shows the value of the downstream Mach number Mdn; the upstream value Mup can be calculated using Equation (11). The finite difference approximations S4 and S5 give good accuracy at sufficiently low Mach numbers, but exhibit significant errors at higher Mach numbers. Once again, the cubic scheme is performs very well along the C characteristic and it is also effective aiong the downstream C+ characteristic at sufficiently small Mach numbers. Unfortunately, these benefits are overshadowed by grossly inadequate behaviour along C+ at high Mach numbers.
Transient polytropic flow
200
803
400
800
1000
Jo0
800
I 1000
-104 200
0.15
0.15WIc+
(vi)c-
0.1
0.1..
3
3 0.05
a
0 I::.::; 200
0.05
400
TE$)
860
1000
3
0 i__?::.m 200 400
TF&
800
1000
Fig. 6. Influence of temperature on the average value of VJkl VAX/D = 0.25, n = 1.4, &,, =0.9; S3, VljV, _ ~~1;S4, Karney and McInnis; S5, cubic; S6, polytropic (exact)].
3.3.3. hJuence of polytropic index. Figure 5 shows the influence of the polytropic index n when all other parameters have the base-case values. With small values of n, heat is received by the fluid as it flows downstream. With larger values, heat is emitted. It is seen that these two extremes have relatively little influence on the accuracy of the finite difference integrations. This result has important implications for extensions of the work to non-polytropic flows. 3.3.4. InfEuence of temperature. Figure 6 shows the influence of the downstream temperature when all other parameters have the base-case values. When the temperatures are small, the wavespeeds are also small and so the velocities implied by any particular Mach number are smaller than those in the base case. The influence of such reductions on the magnitudes of (VI VI),, is partially offset by consequential increases in the times required for waves to travel from one end of the duct to the other. Neither low nor high temperatures have much influence on the accuracy of the approximate methods.
804
A. Vardy and Z. Pan
. Fig. 7. Assumed locations of wavefronts (for integration of environmental terms only).
4. UNSTEADY
POLYTROPIC
FLOW
In transient flows, the axial velocity distribution along an MOC line is not known a priori at arbitrary instants in the numerical solution process. It is not possible for any numerical scheme to obtain exact analytical expressions comparable to those obtained for the steady flow case. Nevertheless, close approximations can be deduced for some special cases and can be used to assess alternative methods of integration. First, all successful schemes must be capable of simulating conditions when only small transients exist anywhere in the pipeline. In the limit, such flows are essentially steady and so the method of integration should reduce to the steady-state case. This is clearly possible with schemes based on Equations (7)-(11). It is not possible at high Mach numbers with approximations such as Equations (13) and (14) or the cubic approximation. Second, successful schemes should be able to reproduce gradually varying conditions resembling uniform acceleration or, more precisely, successive quasi-steady states. Equations (7)-(11) are a good starting point for these flows too. Third, good schemes should be able to model friction in regions of approximately steady flow interspersed with relatively small numbers of strong wavefronts. This type of flow is typified by the classical water-hammer example of a single pipe connecting a water reservoir to a valve. For several pipe periods after sudden valve-closure, the dominant feature is a single wavefront that propagates back and forth along the pipe. Upstream and downstream of this wavefront, the flows are gradually varied, but are approximately steady. In all of these cases, the true skin friction at any particular location depends upon the velocity history as well as upon the instantaneous velocity [9,10]. Nevertheless, attention is focused herein on quasi-steady contributions due to instantaneous velocities only. Additional, inertiadependent contributions may be treated independently when required. 4.1. Idealization In the method described in Section 3 for steady flows, exact integration is possible between end states on characteristic lines because the flow conditions vary continuously in a known manner. In unsteady flows, however, the variation is unknown and there is no simple relationship between the end states. This difficulty must be overcome before use can be made of Equations (7)-(11). In the method of integration adopted herein, each interval is approximated by two regions of steady flow separated by a single wavefront, as illustrated schematically in Fig. 7. On the C+ characteristic, the wavefront is at B and the regions of steady flow are LB= and BAA. Different subscripts are used for these two regions at the point B because this idealization implies that the flow parameters are discontinuous at the wavefront. The idealization is used only for the RHS
Transient polytropic flow
0
0.2
a4
a6
0.8
1
805
0
0.2
0.4
iudn
a15
0
L
0
0.2
1
8 cai
0.1
a a05
0.8
0.15
(iii) c+
3
0.6 Mdn
0.4
Mdn
..
s
a
a6
a8
1
0.06..
01 0
/ a2
a4
0.6
0.8
1
M&l
Fig. 8. Influence of wavefront magnitude on the average value of VIVI [fAx/D = 0.25, n = 1.4, MUP= 0.3. TUP= 300K; S3, VtIV, _ &I; S4, Karney and Mclnnis; S5, cubic; S6, polytropic (exact)].
(environmental) terms in Equation (4), not for the LHS terms, so it does not invalidate the use of the differential Equations (1) and (2). For the purposes of this paper, it is assumed that the locations of the points L and A are known and that the corresponding values of V, p and c, etc., are known. The precise path of the MOC line between the end points cannot (and need not) be known. In a practical implementation, little of this information will be known exactly until the full solution has been obtained at the point A. The overall solution process (location, interpolation, integration) is almost certain to be iterative, especially in high Mach number flows. This is not important herein, however, because attention is focused exclusively on the integration process, not on the method of determining its limit states. The integration of the LHS terms in the MOC equations is described in Section 2.1. When they are of the form shown in Equations (4) and (5), the integration is exact and is independent of idealizations such as that shown in Fig. 7. The proposed integration of the RHS terms is straightforward in principle. The contributions of the two steady-flow regions are determined independently-in the manner described in Section 3-and aire summed to give the required value. There is no contribution from the wavefront because it occupies no space. Mathematically, it is simply a discontinuity between the two regions of integration. The spatial location of the point B is chosen herein at 1/2(xt + xA). Since characteristic lines are rarely straig:ht, the implied temporal location usually differs from 1/2(t,+ tA). Moreover, different times are implied for BA and BL whenever the true flow between A and L is not steady. That is, the implied times are determined by steady-state characteristic lines LB= and BAA.
806
A. Vardy and Z. Pan
Their paths in (x-t) space will differ, in general, from the path in the true unsteady flow. Nevertheless, the total length of duct in which friction is evaluated is exact.
4.2. Integration Unlike the steady-state case, the new scheme should not be regarded as an exact benchmark. Indeed, an infinity of alternative unsteady flow conditions is physically possible between any particular end states and each will have its own value of the integral IV’ dt. In contrast, the values obtained from the numerical methods depend solely upon the end states; the same results are obtained for all intermediate conditions. Nevertheless, it is reasonable to suppose that many unsteady flows will approximate quite closely to the physical model underlying the new method. In these cases, the method may be used as an approximate benchmark. In other cases, no method should be assumed to be accurate. 4.2.1. Choice of 62. In the chosen method of comparison, the same time interval At has been used for all schemes, namely the sum of the intervals obtained for the two regions of integration using polytropic relationships. The polytropic case is evaluated first, giving its result and also giving time intervals such as (rA--rB) and (tB-tL) corresponding to steady polytropic flows in the two regions. The sum of these intervals is then used in the finite difference integrations. This process ensures that differences in Fig. 8 are not attributable to the use of different time steps in the various schemes. In practical applications, there will be differences in the integration time steps and these will influence the values obtained when integrating, say II’” dt. Nevertheless, such differences are not caused by the integration methodology and so they are properly excluded for the purposes of the present paper. In all figures herein, the same length of pipe is used for Cf and C characteristics. This choice, typical of MOC analyses using time-line interpolations, leads to significantly greater time intervals for C lines than for C’ lines. If, instead, equal time intervals were chosen, then different pipe lengths would be implied for the two lines. This would be more typical of analyses using space-line interpolations. 4.2.2. Finite difference uppro.uimutions. Figure 8 compares results obtained using the proposed method with those obtained using the finite difference approximations in Equations (13) and (14). The conditions at the upstream end are held constant (M,,=0.3, Tut,= 300K) and the graphs show the dependence on assumed conditions at the downstream end (0.1 < Mdn < 1.O). No error-graphs are drawn because the exact solution is not known. In the particular case when M,, = Mdn, all four schemes give similar results. Strictly, this case is unsteady because the upstream and downstream Mach numbers would not be equal if the flow were steady. Nevertheless, the Mach number is low and so the most likely real flow is nearly steady. If both M,, and Mdn were equal to some larger value, the most likely real flow would consist of an increasing Mach number from L towards A followed by a reduction and then by a further increase. This condition is implied in the chosen implementation of the polytropic method, the intermediate zone being compressed into a single discontinuity (wavefront) at B.
When the downstream Mach number is smaller than the upstream value (i.e. 0.3 in this example), the finite difference schemes yield smaller values of (VI Vi),” than are obtained with the polytropic approximation. Because the Mach numbers are small, the polytropic idealization approximates to a region 1/2A.u with a Mach number of M,, and a region 1/2Ax with a Mach number of Mdn. The difference between the various results is indicative of different representations of unsteadiness, not compressibility. For example, the polytropic result is typified by 1/2(O.12+O.32) = 0.05 and the V,:V, I A, result is typified by 0.1 x 0.3 = 0.03. When the downstream Mach number is greater than the upstream value, both finite difference schemes again tend to underestimate the most likely values of fV2 dt. Equation (13) yields more probable values than Equation (14) along the C+ characteristic, but less probable ones along C. When Mdn is not much greater than Mupr the implied velocity distributions are broadly similar in all three cases. With larger values of Mdn, however, the finite difference approximations tend to imply smaller velocities in the downstream region than those assumed in the polytropic idealization. The latter is not an exact benchmark, but it may reasonably be inter-
807
Transient polytropic flow
0-b 0
: 0.2
0.4
0.6
0.6
I 1
07
:
0
0.2
:
:
:
0.4
:
0.6
:
:
0.8
1
JfUP
I
(IU)
c+
0.1 ..
3
a 0.05 ..
04 0
I
0.2
:
:
a4
46
:
OS
1
Mup
Fig. 9. Influence of choked flow on the average value of vl Vj VAX/D = 0.25, n = 1.4, Mdn =0.6, T,_+” = 300K; S3, V,lV, _ a,l; S4, Kamey and McInnis; SS, cubic; S6, polytropic (exact)].
preted as giving equal weight to the known conditions at the upstream and downstream ends of the reach. When Mdn exceeds M,,, Karney and Mclnnis’ method (with c = 0.85) gives smaller values than V,IV, + *,I along the downstream characteristic, but larger ones along the upstream characteristic. This is consistent with the steady flow behaviour shown in Figs 3-6, but the implications are different. In Fig. 8, the use oft = 0.85 appears to improve the upstream case, but to worsen the downstream case. In Figs 3-6, the opposite conclusions are obtained. The behaviour aof the cubic approximation is remarkable. Along the C characteristic, it yield results similar to those presented by Kamey and McInnis. Along the C+ characteristic, it agrees closely with the polytropic prediction when Mdn is smaller than about 0.7, but it subsequently predicts decreasing friction with increasing Mdn. This unrealistic behaviour can be traced to the predicted velocity gradient at high Mach numbers. As M --) 1, (aV/ax),+oo and the predicted velocity distribution can include values outside the range from V, to V, + A,. This is one reason why Murray [2] needed to limit the allowable values of the assumed velocity gradients. 4.3. Chokedflow There is one condition that must be treated as a special case. It will rarely occur in computations of real flows, but it can do so when a downstream-moving compression wave creates conditions implying choked flow within the calculation grid length. This case is illustrated in Fig. 9, in which the range of upstream Mach numbers includes values that would not be physically possible except for very short durations. The conditions at the downstream end are held constant (Mdn = 0.6, T,, = 300K).
808
A. Vardy and 2. Pan
When the upstream Mach number is small, the various schemes operate in a broadly similar manner to that shown for the more gentle case considered in Fig. 8, albeit showing greater differences. With successively greater upstream Mach numbers, however, a value is eventually reached for which it is physically impossible for steady flow conditions to prevail for a distance 1/2Ax downstream. In the real flow, the conditions are about to become choked within the integration interval. In Fig. 9, the response to this condition in the polytropic idealization is to evaluate VIVI in the upstream 1/2Ax for an assumed steady flow satisfying the true Mach number at the upstream point L and choked flow at the downstream end A. The skin friction coefficient is reduced to make this possible. To ensure continuous behaviour as the rate of flow increases, the adjustment is made whenever the steady flow asymptote from the point L implies choked conditions upstream of the point A. The finite difference approximations are not affected by the implied choking. Instead, they yield continuously varying results throughout the range considered-even up to an assumed upstream Mach number of 1.0. Physically they are unrealistic, but numerically they are benign. 5. SUMMARY
AND
CONCLUSIONS
A new method of integrating quasi-steady friction in transient polytropic flows has been introduced. It is exact in the particular case of truly steady flows. The method is nearly exact when small disturbances exist in otherwise well-developed steady flows and it is plausible when strong wavefronts exist. The method has been used to assess widely used finite difference methods of integrating friction terms in MOC equations. In steady flows, the approximation VJV, _ a,lAt has been shown to be less accurate than the expression {V, _ At+ 0.85 { Vt-Vl _ al>1V<_ *,lAt} along C+ characteristics, but more accurate along C- characteristics. 4. Because the finite difference integrations are not exact, MOC schemes making use of them will not conserve mass and momentum in asymptotic steady-state conditions. 5. Assuming the new method as an approximate benchmark, the finite difference methods have been shown to underestimate the most likely values of skin friction in regions of significant unsteadiness. When large integration reach& are used, this tendency is- exhibited-even at quite small Mach numbers. 6. In the form presented herein, the new method relies on the availability of exact solutions to steady-state flows. This limits its direct applicability to flows in which only one environmental parameter needs to be integrated.
REFERENCES
1 Karney, B. W. and McInnis, D., Efficient calculation of transient flow in simple pipe networks. Journal of Hydrology Engineeringg, ASCE, 1992, 118(7), 1014-1030. 2. Murray, S. J., Transient analysis of partially full pipe flow. In Proceedings of the International Conference on Unsteady Flow and Fluid Transients, Durham, U.K., 29 September-l October, HR Wallingford. Balkema, Amsterdam, 1992, pp. 143-157. 3. Wylie, E. B., Unsteady internal flows-dimensionless numbers and time constants. In Proceedings of the 7th kernational Conference on Pressure Surges and Fluid Transients in Pipelines and Open Channels, Hairogate, U.K., 16-18 Aoril. BHR GrouD. Mechanical Engineering Publications. 1996. up. 283-288. 4. Rhodes,^F. S., The use of the method of characteristics for the numerical simulation of internal flow discontinuities. Ph.D. thesis, University of Strathclyde, U.K., 1989. 5. Wylie, E. B., Advances in the use of MOC in unsteady pipe flow. Proceedings of the 4th International Conference on Pressure Surges, Bath, U.K., 21-23 September, BHR Group, 1983, pp. 22-37. 6. Vardy, A.E. and Chan, L. I., Rapidly attenuated water-hammer and steel-hammer. In Proceedings of the 4rh International Conference on Pressure Surges, Bath, U.K., 21-23 September, BHR Group, 1983, pp. 22-37. 7. Anderson, A. and Suwan, K., Stability analysis including waterhammer effects. In Proceedings of the International Conference on Unsteady Flow and Fluid Transients, ed. H. R. Wallingford, Durham, U.K., 1992, pp. 71-75. 8. Bergant, A. and Simpson, A. R., Estimating unsteady friction in transient cavitating pipe flow. In Proceedings of the 2nd International Conference on Water Pipeline Systems, Edinburgh, 24-26 May, BHR Group. Mechanical Engineering Publications, 1994, pp. 3-13.
809
Transient polytropic flow
9. Vardy, A. E. and Brown, J. M. B., Transient, turbulent, smooth pipe friction. Journal of Hydrology Research, 1995, 33(4), 435-456.
10. Bughazem, M. B. and Anderson, A., Problems with simple models for damping in unsteady flow. In Proceedings of the 7th International Conference on Pressure Surges and Fluid Transients in Pipelines and Open Channels, Harrogate, U.K., 16-18 April, BHR Group. Mechanical Engineering Publications, 1996, pp. 537-548.
6. APPENDIX 6.1.
Isothermalflow
A (n = 1)
When the flow is isothermal, n = 1 and Equations (4) and (5) are invalid. They may be replaced by
(Al) and When
g=
V-c:
!?!$_~=fg!?+“eu_$~ PA
(0
in which the wavespeed c = J(p/p) is constant. Integration of the LHS terms along the characteristic line LA, for example, leads to
dt =
c{W,) - WL)) + (VA - VL).
(A3)
The integration of the RHS terms is typified along the line LA by Equation (8) which yields A
L
V2&=&
AI-M
2f s L
M
dM (A4)
Similarly, along the line RA,
A
J R
(A5)
7. APPENDIX B 7.1. Evaluation of Mach number at remote location It is frequently necessary to evaluate the Mach number at a specified distance from a flow section at which all conditions are known. In principle, this is achieved by solving Equation (11). In a typical case with a specified value of MA, the only unknown will be the Mach number ML, say. This can be found by solving the equation iteratively. The use of iterative solutions can be time-consuming and this is particularly undesirable when the result is to be used within another iterative process. This can be avoided in the present case by solving the equation in detail at the be inning of each program run, yielding values of (fAx/D) corresponding to a range of Mach numbers from, say, 10-! to 0.995. The values can be stored in arrays for use during the time-marching solution. The method of storing values of (fAx/D) and (M) in the array(s) should be chosen to facilitate their use during the transient solution rather than for convenience of initial evaluation. For the purpose of developing the arrays, however, note that Equation (11) is better suited to deducing distances at which known Mach numbers exist than vice versa.
Computers& Figi& Vol. 26, No. 8, pp. 793-809, 1997 0 1997 ElsevierScience Ltd. All rights reserved Printed in Great Britain 00457930/97 $17.00 t 0.00
PII: s0045-7930(97)00029-7
QUASI-STEADY
FRICTION ALAN
IN TRANSIENT
VARDYt
and ZHIQIU
POLYTROPIC
FLOW
PAN
Civil Engineering Department, University of Dundee, Dundee DDl 4HN, U.K. (Received 24 October 1996; in revised farm 24 October 1997) Abstract-An accurate method of integrating quasi-steady friction in transient polytropic flows is introduced and is justified both physically and analytically. When applied to subsonic flows, it is found to be especially informative at higher Mach numbers where existing methods of integration can lead to serious errors in mass conservation, The method is used herein to assess the accuracy of widely used finite difference representations of quasi-steady friction, attention focusing prin~pally on solutions obtained using the method of characte~sti~. The new method is exact in the special case of steady flows. Of the approximate methods considered, the popular approximation V,/V, _ AJAXis found to perform well. A modification introduced by Karney and McInnis (Journal of Hydology Engineering, ACSE, 1992, 118(7), 1014-1030 [I]) improves the correlation for d’ownstream characteristic lines, but has the opposite effect for upstream characteristics. A simple cubic. scheme, loosely related to a more complex one used successfully by Murray (Proceedings of the International Conference on &steady Flow and Fhdd Transients, 29 September-l October. Balkema, Amsterdam, 1992, pp. 143-157 [2]), illustrates difficulties in higher-order schemes. The new method is not exact in transient Rows, but it gives good accuracy in slowly varying flows and it is plausible in rapidly varying flows. Using it as a benchmark, the approximation k’,lV, _ 4,iAr is again shown to perform well. All of the schemes considered are unreliable when the flow is about to become cho:ked within the region of integration. 0 1997 Elsevier Science Ltd.
NOMENCLATURE A C
C+ c s LHS M MOC n P RHS t V x it Ax Y P f Subscripts A av B dn L, R I r-At up W
cross-sectional area [m*] polytropic wavespeed [=&p/p)] [m s-‘1 characteristic lines in the direction of flow characteristic lines against the direction of flow pipe diameter [m] skin friction efficient [ET,,,/(O.Sp yi Vi)] cross-sectional perimeter [m] left-hand side of MOC equation polytropic Mach number [=V/c] method of characteristics polytropic index pressure [Pa] right-hand side of MOC equation time coordinate [s] average veiocity in cross-section [m P’] distance coordinate [m] elevation [m] time interval (e.g. tA-tL) [s] distance interval (e.g. xA-xL) [m] ratio of principal specific heat capacities mass density kg me31 shear stress [Pa] solution point average value intermediate point on MOC line downstream end of sofution region base points on MOC lines solution time previous calculation time upstream end of solution region wall
~Corresponding author. 793
A. Vardy and Z. Pan
194
1. INTRODUCTION
The method of characteristics (MOC) is the most widely used numerical scheme for simulating transient liquid flows in pipes, especially when there are numerous boundaries. It is capable of yielding accurate predictions for a wide range of flow types provided that limitations on grid sizes are recognized [3]. The picture is less satisfactory for flows in which there are significant variations in wavespeed or where particle velocities are important in comparison with wavespeeds. Common examples include gas flows in long pipelines and free-surface flows in channels. In such cases, the usual methods of expressing MOC equations numerically can lead to serious errors in mass conservation [2,4]. This paper addresses one cause of inaccuracies in MOC when Mach numbers or Froude numbers are significant, namely the integration of terms such as skin friction. One purpose is to introduce a new, highly accurate method of integrating these terms. A more important purpose is to use the new method to assess the validity of simpler, commonly used methods. 1.1. Sources of error in MOC In their simplest form, MOC equations can be integrated exactly. This is so when terms on the right-hand sides of equations such as Equations (4) and (5) herein are zero. It can be possible to obtain exact solutions for simple boundary conditions. In practical applications, the equations must be integrated numerically because of complications from terms describing friction, area change, heat transfer, etc., and from difficulties in locating characteristic lines-especially in fixed-grid methods of solution. Murray [2] set out the three main causes of inaccuracy, namely (i) Location of characteristic lines-the (analytical) properties of the “left-hand side” (LHS) terms in the MOC equations apply only in precisely defined directions in the (x-t) plane. In practice, it is impossible to deduce these directions exactly at all necessary locations. (ii) Interpolation-usually, parameter values must be deduced at locations other than those at which they are notionally known. The resulting interpolation processes involve numerical approximation. (iii) Integration of RHS terms-the numerical integration of the right-hand side (RHS) terms involves approximations when the variations of parameters (e.g. velocity) cannot be expressed analytically. All three of these difficulties can exist in a general method of analysis and all can cause significant errors. Their relative importance is difficult to assess a priori and it is sometimes difficult to distinguish between them. For the purposes of this paper, however, attention is focused exclusively on the third item. Errors due to the first two causes are minimized (often eliminated) by careful selection of the test cases. Where possible, analytical solutions are used to provide exact parameter values at the limits of the region of integration. For simplicity of presentation, the discussion is restricted to polytropic gas flows in horizontal pipes of constant area. In principle, however, the method can be used in many other applications and also in other numerical schemes, e.g. finite differences and finite elements. 2. METHOD
OF CHARACTERISTICS
(MOC)
The equations of continuity and momentum for one-dimensional velocity V along a duct of non-uniform area A may be written as
flow of a perfect gas at a
(1) and
(2) in which p is the mass density, p is the absolute pressure, r, is the wall shear stress, &,, is the length of the duct perimeter, g is the gravitational acceleration and dz/dx is the elevation gradi-
195
Transient polytropic flow
Fig. 1. Characteristic lines
in the (x-r) plane(C’, dx/dt =
V + c; C,
dx/dr = V- c).
ent. Both dz/dx and dA/dx are expressed as an ordinary derivatives to indicate that changes in time are excluded. herein. For the purposes of this paper, all flows are assumed to be polytropic. The effective wave speed c is assumed to satisfy c2 = Y P where n is the polytropic index. The following development is applicable when the flow is not isothermal, that is n # 1. The special case of isothermal flow is presented in Appendix A. Equations (1) and (2) are partial differential equations in two independent variables (x and t). In MOC, they are recast as ordinary differential equations, valid only when specific relationships exist between x and t, namely: dx when -= dt
2
V+c:
dc
dV
(4)
and dx when - = V-c: dt
2 de dV dz --dt=+g;i;;+--7dx. n-ldt
zW&
cVdA (5)
PA
In MOC, the co:nditions dx/dt = I/k c are used to construct characteristic lines in (x-t) space. In Fig. 1, for example, the line LA is chosen so that its gradient at each point satisfies dx/ dt = V + c. As a consequence, Equation (4) is valid at all points on the line and so, in principle, it may be integrated from L to A. Similarly, Equation (5) may be integrated from R to A. 2.1. Integration qf LHS terms The integration of the left-hand side (LHS) terms in Equations (4) and (5) with respect to time t is trivial and the results are analytically exact. Along LA for example:
--$g+g
dt= 1
-& (CA - CL)+ (VA -
VL).
(6)
No assumptions are required in this integration. In particular, no information is required about the manner in w’hich c and V vary along LA (except that they do so continuously). This is the source of the popularity of MOC. It is not always possible to express the equations in a form where the integration can be performed exactly, but it is often possible to do so. In other cases, it m.ight be necessary to use numerical approximations in, for example, coefficients of the term dc/dr.
796
A. Vardy and Z. Pan
0
0.2
0.4
0.6
0.8
1
k/L Fig. 2. Steady polytropic flow in a constant-area duct cfL/O = 50, n = 1.4, pdn = 100 kPa, Tdn= 300K).
2.2. Integration of RHS terms
The integration of the RHS terms in Equations (4) and (5) is intrinsically more difficult than that of the LHS terms because of the existence of expressions such as cV and VIVI. To integrate these exactly would require a knowledge of their variations in time, which is not possible a priori. In practice, therefore, the integration must be undertaken numerically. There are two circumstances in which the need for numerical evaluation does not have a serious adverse influence on the final solution. The first is when the magnitudes of the RHS terms are much smaller than the magnitudes of the LHS terms. This is quite common in highly transient, low Mach number flows. It is less common in higher Mach number flows, partly because transient events disperse more readily and partly because velocity-dependent terms tend to be much larger. The second circumstance is when the true variation of, say, VIVJ is well represented by the numerical approximation. This is quite common in regions of low Mach number flow remote from transient disturbances. In such regions, there is little variation in velocity on any particular characteristic line and so all reasonable approximations to VIVI yield similar values. At higher Mach numbers, however, considerable variations in velocity can exist, even in a steady flow, so different methods of approximating VIVI can yield significantly different values. One purpose of this paper is to extend the range of cases when the numerical integration can be performed accurately-by removing the restriction to low Mach number flows. This is achieved by providing a good approximation to the true variation of velocity along characteristic lines. 2.3. Qualitative d@erence between LHS and RHS terms The distinction between the LHS and RHS terms in Equations (4) and (5) is emphasized herein by referring to the latter as environmental terms. These describe background conditions that exist even in the absence of waves (i.e. in a steady flow). Mathematically, they may be functions of the dependent parameters, but they must not depend on derivatives of these parameters. Physically, they are more closely related to steady-flow considerations than to unsteady ones. Small waves tend to have little influence on them. In contrast, the LHS terms are strongly dependent on unsteadiness; they can be large even when the RHS terms are zero. 3. STEADY
POLYTROPIC
FLOW
The detailed development begins with the special case of steady flow. This is important in its own right because many unsteady flows ultimately become steady. It is also the condition in
Transient polytropic flow
797
which lack of m,ass conservation is most easily detected. The steady flow case has even greater significance herem, however, because it underpins the methodology adopted for unsteady flows in Section 4.
3.1. Integration between prescribed Mach numbers Figure 2 shows the axial velocity dist~bution in a steady, polytropic flow of air along a horizontal duct of constant area and constant friction coefficient f (~~~~~0.5~~1~1~). Near the upstream end (low Mach numbers), the velocity increases gradually as the reduction in pressure due to friction causes a decrease in density. Near the downstream end (high Mach numbers), the rates of change are much greater because the skin friction stresses are greater. At any position along the duct, the local rate of change of the polytropic Mach number M with distance x satisfies dM S(n+l) dx=D
N3 (l-M*)’
(7)
Using this expression, the integral s,” I/ * dt may be developed along a characteristic line satisfying dx/dt = Y i- c as
JA~x
=
,I_ l+M
dx
R 1-M
D
=f(n
J LCM
(8)
d”.
To integrate the last of these expressions exactly, it is necessary to eliminate the wavespeed c. The method used to do this depends upon the value of the polytropic index n. The particular case of isothermal flow (n = 1) is considered in Appendix A. The procedure for other values of n is as follows. For polytrop.ic flows, the product c@ - ‘)‘w + ‘) is constant. This enables the wavespeed c inside the integral in Equation (7) to be eliminated, giving n-l A
J
ff*&
L
n+l DcAM~
A
1)
I_L
f(n+
=
I-M
2n Mn+l
d&I
Similarly, along the characteristic line RA n-l
J
A
R
V*dt=-
DC&~ f@+l)
A
1+&f
J R --3TdM _
Mn+l
A. Vardy and Z. Pan
798
=-
DQM;+~ f
3.2. Integration between prescribed locations
The preceding analysis enables s;’ V 2 dt and f,” V2 dt to be evaluated between flow sections at which the Mach numbers are known. In practical applications, however, it is more usual to integrate between prescribed locations than between prescribed Mach numbers. In such cases, the Mach number is known at one end of the region, but it is not known a priori at the other end. Fortunately, the determination of the other Mach number is straightforward. Its value is obtained by integrating Equation (7), giving
The existence of the logarithm in this expression complicates the numerical determination of the unknown Mach number because it necessitates an iterative method of solution. This is undesirable for present purposes because the integration is performed many times in the proposed method of analysis for transient flows (see Section 4). A simple way of overcoming this complication is outlined in Appendix B. 3.3. Comparison with finite d@erence approximations Equations (7)-(11) are exact. Inaccuracies can arise when they are used in MOC analysesbecause of inaccurate values of iwL or xA-xL, say (see Section l.l)-but the integrations themselves are exact. They may therefore be used to assess the accuracy of commonly used approximations. The most common ways of integrating environmental terms utilize finite difference approximations such as (along LA): A JL
V] VI dtx( Vi V&“At.
02)
in which At = (tA--tL) and, by implication, the integration is carried out along a characteristic line. Many ways of approximating (VI VI),, are possible. Examples include (V, _ AtlVr _ J), (V,l Vg/) and (V,,l V,,i) where Vay= 1/2(V, _ dt+ k;). The choices are essentially arbitrary. It cannot be known a priori which will give the best result. Two approximations used extensively herein are (VlVl),,AtX
k’,Iv,-AttAt
(13)
and (VI VI)&=
f V;-Ar
+ @“-t -
vi-At))
f v,-AdAt.
(14)
The approximation VJ V, _ ~~1was presented by Wylie [5] and by Vardy and Chan [6] at the same conference. It is simple, easy to code and it does not necessitate an iterative method of solution. Most importantly, it appears to be unconditionally stable [5,7,8]. Equation (14) was introduced by Karney and McInnis [l]. It is a weighted sum of Equation (13) and the simple approximation V, _ &I Vt _ &it. When c = 1.0, it becomes identical to Equation (13). Karney and McInnis tentatively recommended the use of t; = 0.85, This value is implied herein whenever Equation (14) is used. Murray [2] recognized that the integration could be performed more accurately if values of the velocity could be estimated at intermediate positions on the characteristic lines. The inte-
Transient polytropic flow
0.15
I
(v)
799
I
c+
0.1
0.1
3 il
0 a 0.05
0.05
1
,t=L 0
0.1
0
0.2
0.3
0.4
0.5
f_ Fig. 3. Influence of fAx/D on the average value of Y(vl [n = 1.4, Mdn= 0.9, Tdn= 300K; Sl, v, _ & v, _ A,[; S2, v,l v,l; S3, v,l v, _ ~~1;S4, Karney and McInnis; SS, cubic; S6, polytropic (exact)].
gration could then be carried out in, say, four intervals of 1/4At instead of one interval of At. To obtain the intermediate velocities, he first located the x-positions of the 1/4At points, using dx/dt = Vf c at the ends L and A to obtain a cubic x(t). He then deduced the corresponding depths (JJ) in his free-surface flow, using estimated values of dy/dt at the ends L and A to obtain a cubic y(x). Finally, the corresponding velocities were obtained, based on an assumed linear distribution of fIowrate. Murray’s method is not applied herein. In part, this is because his description intentionally left some decisions open (he wanted to propose a methodology rather than define a universally applicable imple:mentation). In addition, it is because the method involves adjustments that give rise to discontinuous behaviour when presented in the manner used herein. Instead of applying Murray’s method directly, a cubic polynomial V(t) is created using (a) V and t at the ends of the characteristic lines and (b) assumed velocity gradients iilL’/& at these locations. The assumed gradients satisfy
A. Vardy and Z. Pan
800
Table 1. Base values for comparisons of CV*dr 0.25 1.4 0.9 300K
Skin friction fAx/D Polytropic index n Downstream Mach number Mdn Downstream temperature Tdn
at=[ax 1 s dt av
av
--
dx
(15)
where [aV/ax], denotes the gradient that would exist in a steady flow with the same velocity and dx/dt is the gradient of the appropriate characteristic line. Murray began with these gradients, but did not use [aV/ax], directly. Instead, he first made a plausible adjustment and then limited the allowable values to avoid (a) inflexions and (b) the implied existence, within the range t - At to t, of velocities outwith the range V! _ At to V,. These adjustments are not made herein, thus simplifying the presentation, but at the expense of failing to constrain undesirable behaviour when the Mach number approaches unity, implying [aV/ax],+oo. Having obtained a cubic polynomial satisfying the velocities and the gradients [aV/ax],, the evaluation of the integral JV’ dt between the end states at, say, L and A is straightforward. In the remainder of this paper, this is referred to simply as the ‘cubic’ scheme. 3.3.1. Influence offsx/d. Figure 3 illustrates the influence of the parameter fAx/D on the accuracy of five methods of approximation (defined in the figure caption). The top two boxes show the average values of VI VI implied by the various schemes along the characteristic lines LA and RA, respectively. The required values of the integrals $ V 2 dt and s,” V 2 dt are the products of (VI VI),, and the time intervals At z (t,-tl _ A, ) shown in the bottom two boxes. The interval tA - tR consistently exceeds tA-tL because the simulated direction of flow is always from L to R. Upstream-moving waves travel more slowly than downstream-moving waves. The middle pair of boxes shows percentage errors associated with three of these schemes. In all figures herein, the polytropic case is shown as a continuous line. In the lower boxes, it coincides with error = 0%. Except for fAx/D, all parameters in Fig. 3 have the base-case values listed in Table 1 (e.g. n = 1.4). Values for the polytropic (exact) case are obtained for the C+ characteristic line LA as follows: (i) given fAx/D and the values of all parameters at the (downstream) end A, evaluate the interval At = tA-tL by integrating dt = dx/( V i- c); (ii) evaluate all parameters at the upstream end L using Equation (11) and polytropic state relationships; (iii) evaluate jf V 2 dt using Equation (9); (iv) divide by At to give (VI Vl)aY. A similar procedure is used for the C characteristic line RA, for which the point A is chosen at the upstream end of the same duct. All results shown in Fig. 3 have been obtained using exact values of all parameters at the various limits of integration. That is, all end states satisfy analytical relationships and the paths of the characteristic lines are exact. This ensures a fair comparison of the various methods of integration, with no complications arising from location or interpolation errors. Schemes Sl and S2 use values of II VI determined at only one end of the duct. They are seen to be significantly less accurate than the other combinations and so they are not considered hereafter. Because of their relevance to Scheme S5 [Equation (14)], however, note that each overestimates along one characteristic, but underestimates along the other. The curves obtained using Equations (13) and (14) (with E = 0.85) are the most accurate of the approximations considered in the figure, but the errors are significant when fAx/D is sufficiently large. Expressed as a percentage of (VI VI),, , however, the error resulting from the use of Equation (13) is smaller along the C characteristic than along the C+ characteristic. The absolute errors in the integral JV” dt are approximately equal along the two characteristics, how-
Transient polytropic flow
0
0.2
0.4
0.6
0.8
801
0
1
0.2
0
0.2
0.4
0.6
0.8
0.4
0.6
0.8
1
M
M
1
-10I 0
0.2
M
0.4
0.6
0.8
1
M
0.15 8
0
0.2
0.4
0.6 M
0.8
1
0
c-
0.2
I
0.4
0.6
0.8
1
M
Fig. 4. Influence of Mach number on the average value of v V’j[fAx/D = 0.25, n = 1.4, Td,,= 300K; S3, V,l V, _ *,I; S4, Karney and McInnis; S5, cubic; S6, polytropic (exact)].
ever, because of the influence of the time interval At. In contrast, Equation (14) gives approximately equal percentage errors along the two characteristics, but the absolute error jV* dt is much larger along the C characteristic than along the C+ characteristic. Thus the improvement reported by Karney and McInnis for transient, low Mach number flows is not obtained in steady, higher Mach number flows. For all steady flows considered herein, it would be useful to choose 6 < 1 along downstream characteristics, but t: > 1 along upstream characteristics. The detailed quantification of this process would merit a full investigation in its own right because the application of this elegant method is so simple. The cubic scheme performs poorly on the downstream characteristic, but very well on the upstream characteristic. These different behaviours are attributable to the different influences of the gradients of the characteristic lines used in Equation (1 S), but the authors have been unable to find a satisfactory way of exploiting this fact. Adjustments such as those proposed by
802
A. Vardy and 2. Pan
1
1.2
1.4
1.6
1.8
2
1.4
1.2
1
30 1
I
1.6
1.8
R
II
30 ,
2o._.. “‘...I”...... ..^ .. ..YL . l+.._._ _........ I
I
c-
I
il______~ ____.~
;_!
-20
(iv)
2
-- -- --,__,__,______\__-_ *- -- ,______.__.__,__.__
3.
1.2
1
1.4
1.6
1.8
2
II
1
1.2
1.4
1.6
1.8
II
0.15 /
1
2
0.15 (vi) c+-
0.1 -. $ i5
0.M -.
01 1
:
: 1.2
:
: 1.4
: n
: 1.6
: 1.8
:
1 2
04 1
:
1.2
:
:
:
1.4
:
:
1.6
;
1.8
2
n
Fig. 5. Influence of the polytropic index on the average value of fl vi @&.x/B= 0.25, k&=0.9, T& = 300K; 53, YtlV, _ &I; 54, Karney and McInnis; S.5, cubic; 56, polytropic (exact)].
Murray [2] have been found to reduce the errors on the downstream characteristic, but their application renders the process cumbersome and the integral does not necessarily increase monotonically with ~&X/D. Moreover, in the authors’ implementations, the same adjustments had a detrimental influence on the integration along the upstream characteristic. 3.3.2. Injhence of Much number. Figure 4 shows the influence of the Mach number when all other parameters have the base-case values listed in Table 1. The abscissa shows the value of the downstream Mach number Mdn; the upstream value Mup can be calculated using Equation (11). The finite difference approximations S4 and S5 give good accuracy at sufficiently low Mach numbers, but exhibit significant errors at higher Mach numbers. Once again, the cubic scheme is performs very well along the C characteristic and it is also effective aiong the downstream C+ characteristic at sufficiently small Mach numbers. Unfortunately, these benefits are overshadowed by grossly inadequate behaviour along C+ at high Mach numbers.
Transient polytropic flow
200
803
400
800
1000
Jo0
800
I 1000
-104 200
0.15
0.15WIc+
(vi)c-
0.1
0.1..
3
3 0.05
a
0 I::.::; 200
0.05
400
TE$)
860
1000
3
0 i__?::.m 200 400
TF&
800
1000
Fig. 6. Influence of temperature on the average value of VJkl VAX/D = 0.25, n = 1.4, &,, =0.9; S3, VljV, _ ~~1;S4, Karney and McInnis; S5, cubic; S6, polytropic (exact)].
3.3.3. hJuence of polytropic index. Figure 5 shows the influence of the polytropic index n when all other parameters have the base-case values. With small values of n, heat is received by the fluid as it flows downstream. With larger values, heat is emitted. It is seen that these two extremes have relatively little influence on the accuracy of the finite difference integrations. This result has important implications for extensions of the work to non-polytropic flows. 3.3.4. InfEuence of temperature. Figure 6 shows the influence of the downstream temperature when all other parameters have the base-case values. When the temperatures are small, the wavespeeds are also small and so the velocities implied by any particular Mach number are smaller than those in the base case. The influence of such reductions on the magnitudes of (VI VI),, is partially offset by consequential increases in the times required for waves to travel from one end of the duct to the other. Neither low nor high temperatures have much influence on the accuracy of the approximate methods.
804
A. Vardy and Z. Pan
. Fig. 7. Assumed locations of wavefronts (for integration of environmental terms only).
4. UNSTEADY
POLYTROPIC
FLOW
In transient flows, the axial velocity distribution along an MOC line is not known a priori at arbitrary instants in the numerical solution process. It is not possible for any numerical scheme to obtain exact analytical expressions comparable to those obtained for the steady flow case. Nevertheless, close approximations can be deduced for some special cases and can be used to assess alternative methods of integration. First, all successful schemes must be capable of simulating conditions when only small transients exist anywhere in the pipeline. In the limit, such flows are essentially steady and so the method of integration should reduce to the steady-state case. This is clearly possible with schemes based on Equations (7)-(11). It is not possible at high Mach numbers with approximations such as Equations (13) and (14) or the cubic approximation. Second, successful schemes should be able to reproduce gradually varying conditions resembling uniform acceleration or, more precisely, successive quasi-steady states. Equations (7)-(11) are a good starting point for these flows too. Third, good schemes should be able to model friction in regions of approximately steady flow interspersed with relatively small numbers of strong wavefronts. This type of flow is typified by the classical water-hammer example of a single pipe connecting a water reservoir to a valve. For several pipe periods after sudden valve-closure, the dominant feature is a single wavefront that propagates back and forth along the pipe. Upstream and downstream of this wavefront, the flows are gradually varied, but are approximately steady. In all of these cases, the true skin friction at any particular location depends upon the velocity history as well as upon the instantaneous velocity [9,10]. Nevertheless, attention is focused herein on quasi-steady contributions due to instantaneous velocities only. Additional, inertiadependent contributions may be treated independently when required. 4.1. Idealization In the method described in Section 3 for steady flows, exact integration is possible between end states on characteristic lines because the flow conditions vary continuously in a known manner. In unsteady flows, however, the variation is unknown and there is no simple relationship between the end states. This difficulty must be overcome before use can be made of Equations (7)-(11). In the method of integration adopted herein, each interval is approximated by two regions of steady flow separated by a single wavefront, as illustrated schematically in Fig. 7. On the C+ characteristic, the wavefront is at B and the regions of steady flow are LB= and BAA. Different subscripts are used for these two regions at the point B because this idealization implies that the flow parameters are discontinuous at the wavefront. The idealization is used only for the RHS
Transient polytropic flow
0
0.2
a4
a6
0.8
1
805
0
0.2
0.4
iudn
a15
0
L
0
0.2
1
8 cai
0.1
a a05
0.8
0.15
(iii) c+
3
0.6 Mdn
0.4
Mdn
..
s
a
a6
a8
1
0.06..
01 0
/ a2
a4
0.6
0.8
1
M&l
Fig. 8. Influence of wavefront magnitude on the average value of VIVI [fAx/D = 0.25, n = 1.4, MUP= 0.3. TUP= 300K; S3, VtIV, _ &I; S4, Karney and Mclnnis; S5, cubic; S6, polytropic (exact)].
(environmental) terms in Equation (4), not for the LHS terms, so it does not invalidate the use of the differential Equations (1) and (2). For the purposes of this paper, it is assumed that the locations of the points L and A are known and that the corresponding values of V, p and c, etc., are known. The precise path of the MOC line between the end points cannot (and need not) be known. In a practical implementation, little of this information will be known exactly until the full solution has been obtained at the point A. The overall solution process (location, interpolation, integration) is almost certain to be iterative, especially in high Mach number flows. This is not important herein, however, because attention is focused exclusively on the integration process, not on the method of determining its limit states. The integration of the LHS terms in the MOC equations is described in Section 2.1. When they are of the form shown in Equations (4) and (5), the integration is exact and is independent of idealizations such as that shown in Fig. 7. The proposed integration of the RHS terms is straightforward in principle. The contributions of the two steady-flow regions are determined independently-in the manner described in Section 3-and aire summed to give the required value. There is no contribution from the wavefront because it occupies no space. Mathematically, it is simply a discontinuity between the two regions of integration. The spatial location of the point B is chosen herein at 1/2(xt + xA). Since characteristic lines are rarely straig:ht, the implied temporal location usually differs from 1/2(t,+ tA). Moreover, different times are implied for BA and BL whenever the true flow between A and L is not steady. That is, the implied times are determined by steady-state characteristic lines LB= and BAA.
806
A. Vardy and Z. Pan
Their paths in (x-t) space will differ, in general, from the path in the true unsteady flow. Nevertheless, the total length of duct in which friction is evaluated is exact.
4.2. Integration Unlike the steady-state case, the new scheme should not be regarded as an exact benchmark. Indeed, an infinity of alternative unsteady flow conditions is physically possible between any particular end states and each will have its own value of the integral IV’ dt. In contrast, the values obtained from the numerical methods depend solely upon the end states; the same results are obtained for all intermediate conditions. Nevertheless, it is reasonable to suppose that many unsteady flows will approximate quite closely to the physical model underlying the new method. In these cases, the method may be used as an approximate benchmark. In other cases, no method should be assumed to be accurate. 4.2.1. Choice of 62. In the chosen method of comparison, the same time interval At has been used for all schemes, namely the sum of the intervals obtained for the two regions of integration using polytropic relationships. The polytropic case is evaluated first, giving its result and also giving time intervals such as (rA--rB) and (tB-tL) corresponding to steady polytropic flows in the two regions. The sum of these intervals is then used in the finite difference integrations. This process ensures that differences in Fig. 8 are not attributable to the use of different time steps in the various schemes. In practical applications, there will be differences in the integration time steps and these will influence the values obtained when integrating, say II’” dt. Nevertheless, such differences are not caused by the integration methodology and so they are properly excluded for the purposes of the present paper. In all figures herein, the same length of pipe is used for Cf and C characteristics. This choice, typical of MOC analyses using time-line interpolations, leads to significantly greater time intervals for C lines than for C’ lines. If, instead, equal time intervals were chosen, then different pipe lengths would be implied for the two lines. This would be more typical of analyses using space-line interpolations. 4.2.2. Finite difference uppro.uimutions. Figure 8 compares results obtained using the proposed method with those obtained using the finite difference approximations in Equations (13) and (14). The conditions at the upstream end are held constant (M,,=0.3, Tut,= 300K) and the graphs show the dependence on assumed conditions at the downstream end (0.1 < Mdn < 1.O). No error-graphs are drawn because the exact solution is not known. In the particular case when M,, = Mdn, all four schemes give similar results. Strictly, this case is unsteady because the upstream and downstream Mach numbers would not be equal if the flow were steady. Nevertheless, the Mach number is low and so the most likely real flow is nearly steady. If both M,, and Mdn were equal to some larger value, the most likely real flow would consist of an increasing Mach number from L towards A followed by a reduction and then by a further increase. This condition is implied in the chosen implementation of the polytropic method, the intermediate zone being compressed into a single discontinuity (wavefront) at B.
When the downstream Mach number is smaller than the upstream value (i.e. 0.3 in this example), the finite difference schemes yield smaller values of (VI Vi),” than are obtained with the polytropic approximation. Because the Mach numbers are small, the polytropic idealization approximates to a region 1/2A.u with a Mach number of M,, and a region 1/2Ax with a Mach number of Mdn. The difference between the various results is indicative of different representations of unsteadiness, not compressibility. For example, the polytropic result is typified by 1/2(O.12+O.32) = 0.05 and the V,:V, I A, result is typified by 0.1 x 0.3 = 0.03. When the downstream Mach number is greater than the upstream value, both finite difference schemes again tend to underestimate the most likely values of fV2 dt. Equation (13) yields more probable values than Equation (14) along the C+ characteristic, but less probable ones along C. When Mdn is not much greater than Mupr the implied velocity distributions are broadly similar in all three cases. With larger values of Mdn, however, the finite difference approximations tend to imply smaller velocities in the downstream region than those assumed in the polytropic idealization. The latter is not an exact benchmark, but it may reasonably be inter-
807
Transient polytropic flow
0-b 0
: 0.2
0.4
0.6
0.6
I 1
07
:
0
0.2
:
:
:
0.4
:
0.6
:
:
0.8
1
JfUP
I
(IU)
c+
0.1 ..
3
a 0.05 ..
04 0
I
0.2
:
:
a4
46
:
OS
1
Mup
Fig. 9. Influence of choked flow on the average value of vl Vj VAX/D = 0.25, n = 1.4, Mdn =0.6, T,_+” = 300K; S3, V,lV, _ a,l; S4, Kamey and McInnis; SS, cubic; S6, polytropic (exact)].
preted as giving equal weight to the known conditions at the upstream and downstream ends of the reach. When Mdn exceeds M,,, Karney and Mclnnis’ method (with c = 0.85) gives smaller values than V,IV, + *,I along the downstream characteristic, but larger ones along the upstream characteristic. This is consistent with the steady flow behaviour shown in Figs 3-6, but the implications are different. In Fig. 8, the use oft = 0.85 appears to improve the upstream case, but to worsen the downstream case. In Figs 3-6, the opposite conclusions are obtained. The behaviour aof the cubic approximation is remarkable. Along the C characteristic, it yield results similar to those presented by Kamey and McInnis. Along the C+ characteristic, it agrees closely with the polytropic prediction when Mdn is smaller than about 0.7, but it subsequently predicts decreasing friction with increasing Mdn. This unrealistic behaviour can be traced to the predicted velocity gradient at high Mach numbers. As M --) 1, (aV/ax),+oo and the predicted velocity distribution can include values outside the range from V, to V, + A,. This is one reason why Murray [2] needed to limit the allowable values of the assumed velocity gradients. 4.3. Chokedflow There is one condition that must be treated as a special case. It will rarely occur in computations of real flows, but it can do so when a downstream-moving compression wave creates conditions implying choked flow within the calculation grid length. This case is illustrated in Fig. 9, in which the range of upstream Mach numbers includes values that would not be physically possible except for very short durations. The conditions at the downstream end are held constant (Mdn = 0.6, T,, = 300K).
808
A. Vardy and 2. Pan
When the upstream Mach number is small, the various schemes operate in a broadly similar manner to that shown for the more gentle case considered in Fig. 8, albeit showing greater differences. With successively greater upstream Mach numbers, however, a value is eventually reached for which it is physically impossible for steady flow conditions to prevail for a distance 1/2Ax downstream. In the real flow, the conditions are about to become choked within the integration interval. In Fig. 9, the response to this condition in the polytropic idealization is to evaluate VIVI in the upstream 1/2Ax for an assumed steady flow satisfying the true Mach number at the upstream point L and choked flow at the downstream end A. The skin friction coefficient is reduced to make this possible. To ensure continuous behaviour as the rate of flow increases, the adjustment is made whenever the steady flow asymptote from the point L implies choked conditions upstream of the point A. The finite difference approximations are not affected by the implied choking. Instead, they yield continuously varying results throughout the range considered-even up to an assumed upstream Mach number of 1.0. Physically they are unrealistic, but numerically they are benign. 5. SUMMARY
AND
CONCLUSIONS
A new method of integrating quasi-steady friction in transient polytropic flows has been introduced. It is exact in the particular case of truly steady flows. The method is nearly exact when small disturbances exist in otherwise well-developed steady flows and it is plausible when strong wavefronts exist. The method has been used to assess widely used finite difference methods of integrating friction terms in MOC equations. In steady flows, the approximation VJV, _ a,lAt has been shown to be less accurate than the expression {V, _ At+ 0.85 { Vt-Vl _ al>1V<_ *,lAt} along C+ characteristics, but more accurate along C- characteristics. 4. Because the finite difference integrations are not exact, MOC schemes making use of them will not conserve mass and momentum in asymptotic steady-state conditions. 5. Assuming the new method as an approximate benchmark, the finite difference methods have been shown to underestimate the most likely values of skin friction in regions of significant unsteadiness. When large integration reach& are used, this tendency is- exhibited-even at quite small Mach numbers. 6. In the form presented herein, the new method relies on the availability of exact solutions to steady-state flows. This limits its direct applicability to flows in which only one environmental parameter needs to be integrated.
REFERENCES
1 Karney, B. W. and McInnis, D., Efficient calculation of transient flow in simple pipe networks. Journal of Hydrology Engineeringg, ASCE, 1992, 118(7), 1014-1030. 2. Murray, S. J., Transient analysis of partially full pipe flow. In Proceedings of the International Conference on Unsteady Flow and Fluid Transients, Durham, U.K., 29 September-l October, HR Wallingford. Balkema, Amsterdam, 1992, pp. 143-157. 3. Wylie, E. B., Unsteady internal flows-dimensionless numbers and time constants. In Proceedings of the 7th kernational Conference on Pressure Surges and Fluid Transients in Pipelines and Open Channels, Hairogate, U.K., 16-18 Aoril. BHR GrouD. Mechanical Engineering Publications. 1996. up. 283-288. 4. Rhodes,^F. S., The use of the method of characteristics for the numerical simulation of internal flow discontinuities. Ph.D. thesis, University of Strathclyde, U.K., 1989. 5. Wylie, E. B., Advances in the use of MOC in unsteady pipe flow. Proceedings of the 4th International Conference on Pressure Surges, Bath, U.K., 21-23 September, BHR Group, 1983, pp. 22-37. 6. Vardy, A.E. and Chan, L. I., Rapidly attenuated water-hammer and steel-hammer. In Proceedings of the 4rh International Conference on Pressure Surges, Bath, U.K., 21-23 September, BHR Group, 1983, pp. 22-37. 7. Anderson, A. and Suwan, K., Stability analysis including waterhammer effects. In Proceedings of the International Conference on Unsteady Flow and Fluid Transients, ed. H. R. Wallingford, Durham, U.K., 1992, pp. 71-75. 8. Bergant, A. and Simpson, A. R., Estimating unsteady friction in transient cavitating pipe flow. In Proceedings of the 2nd International Conference on Water Pipeline Systems, Edinburgh, 24-26 May, BHR Group. Mechanical Engineering Publications, 1994, pp. 3-13.
809
Transient polytropic flow
9. Vardy, A. E. and Brown, J. M. B., Transient, turbulent, smooth pipe friction. Journal of Hydrology Research, 1995, 33(4), 435-456.
10. Bughazem, M. B. and Anderson, A., Problems with simple models for damping in unsteady flow. In Proceedings of the 7th International Conference on Pressure Surges and Fluid Transients in Pipelines and Open Channels, Harrogate, U.K., 16-18 April, BHR Group. Mechanical Engineering Publications, 1996, pp. 537-548.
6. APPENDIX 6.1.
Isothermalflow
A (n = 1)
When the flow is isothermal, n = 1 and Equations (4) and (5) are invalid. They may be replaced by
(Al) and When
g=
V-c:
!?!$_~=fg!?+“eu_$~ PA
(0
in which the wavespeed c = J(p/p) is constant. Integration of the LHS terms along the characteristic line LA, for example, leads to
dt =
c{W,) - WL)) + (VA - VL).
(A3)
The integration of the RHS terms is typified along the line LA by Equation (8) which yields A
L
V2&=&
AI-M
2f s L
M
dM (A4)
Similarly, along the line RA,
A
J R
(A5)
7. APPENDIX B 7.1. Evaluation of Mach number at remote location It is frequently necessary to evaluate the Mach number at a specified distance from a flow section at which all conditions are known. In principle, this is achieved by solving Equation (11). In a typical case with a specified value of MA, the only unknown will be the Mach number ML, say. This can be found by solving the equation iteratively. The use of iterative solutions can be time-consuming and this is particularly undesirable when the result is to be used within another iterative process. This can be avoided in the present case by solving the equation in detail at the be inning of each program run, yielding values of (fAx/D) corresponding to a range of Mach numbers from, say, 10-! to 0.995. The values can be stored in arrays for use during the time-marching solution. The method of storing values of (fAx/D) and (M) in the array(s) should be chosen to facilitate their use during the transient solution rather than for convenience of initial evaluation. For the purpose of developing the arrays, however, note that Equation (11) is better suited to deducing distances at which known Mach numbers exist than vice versa.