Vardy-Tijsseling 2015
Vardy-Tijsseling 2015
Vardy-Tijsseling 2015
net/publication/292980982
CITATIONS READS
7 451
2 authors, including:
Arris S Tijsseling
Eindhoven University of Technology
121 PUBLICATIONS 1,766 CITATIONS
SEE PROFILE
All content following this page was uploaded by Arris S Tijsseling on 04 February 2016.
A E Vardy
School of Science & Engineering, University of Dundee, UK
A S Tijsseling
Department of Mathematics and Computer Science
Eindhoven University of Technology, The Netherlands
ABSTRACT
The method of characteristics (MOC) has undoubtedly been a wonderful bonus for the
pressure surge community. Its simplicity is almost too good to believe and, when used
properly, the accuracy of its predictions in practical applications is a great comfort.
When the method is applied to liquid flows in elastic pipes, these attributes are so robust
that they are commonly taken for granted. This forum paper explores how far we may
stray from “liquid flows in elastic pipes” before reliance on the attributes ceases to be
safe. This is done by independent assessments of (i) LHS terms – i.e. differential terms,
(ii) RHS terms – e.g. friction, mass addition and (iii) the gradients of the characteristics.
Questions that are answered include:
• Is it possible to express MOC equations with truly constant coefficients of LHS terms?
• Is it possible to show how quasi-steady friction influences wave speed?
• Is it possible to be confident of reproducing steady flow correctly?
NOMENCLATURE
c sound speed
CRHS RHS terms in the continuity equation
g gravitational acceleration
K effective bulk modulus
LHS left hand side
MRHS RHS terms in the continuity equation
MOC method of characteristics
p pressure
RHS right hand side
t time coordinate
U mean velocity
x distance coordinate
Greek symbols
α coefficients of RHS terms in compatibility equations
β coefficients of RHS terms in compatibility equations
Δ finite interval
ρ fluid density
The paper is intended to promote – or even provoke – discussion. The authors hope that
it does not contain serious errors, but acknowledge that it pushes boundaries a little and
so emphasize that it should not be regarded as an attempt to provide formal methods for
widespread use. We hope that it will not be used for that purpose, but will be pleased if it
merits future citations as a contributor to more rigorous follow-up work. All readers will
already have a good appreciation of much of the paper’s content. We suspect however,
that few will have full understanding of every issue raised. Indeed, the authors do not
themselves claim to be members of such a rare breed. What we bring to the party is a
slightly unusual background in MOC. The first author has spent most of his professional
career working on unsteady gas flows and his best known academic contribution is in
unsteady skin friction. The second author has had a longstanding interest in fluid-
structure interaction (FSI) where waves propagate at greatly different speeds in the fluid
and solid. Also, he is as familiar with the analysis of unsteady flows in the frequency
domain as he is with MOC’s time domain approach.
The main body of the paper begins with a presentation of the MOC equations that will be
familiar to most readers. The particular form of the presentation is straightforward, but it
provides a basis for the remainder of the paper. In particular, it avoids a need for
annoying repetition in subsequent sections. There are then three sections discussing
aspects of:
• The left hand sides (LHS) of the compatibility equations – i.e. the differential terms
• The right hand side (RHS) of the compatibility equations – i.e. the non-differential
terms
• The characteristic directions
The paper then has some closing remarks. It would be presumptuous to claim these as
formal conclusions so we do not do so. It was our original intention to have a fourth core
section dealing with boundary conditions. However, these have such a profound
influence on MOC outcomes that they merit fuller attention in their own right. Also,
their dominant importance would detract from the somewhat subtler points that are the
focus herein. Accordingly, that work has been excised and will be addressed in a
different manner, perhaps using other opportunities than conferences present.
The whole of the paper is written in a somewhat informal style that contrasts with the
authors’ usual papers. This reflects its declared purpose of provoking discussion. It also
helps to protect readers from assuming that the paper tells a complete story. It doesn’t; it
simply shows where careful thought is needed by anyone who wishes to stray outside the
boundaries where conventional water hammer is so hugely successful.
(2.1)
+ + =
in which p, ρ & U denote the pressure, mean velocity & density and x,t are the spatial and
time co-ordinates. The terms CRHS and MRHS depend upon the nature of the particular
case under study. For instance, CRHS can allow for spatially varied flow and MRHS can
additionally allow for friction and for body forces (e.g. gravitational). If, as usual in
water-hammer, these equations are simplified by neglecting the convective terms U∂/∂x
in comparison with inertial terms ∂/∂t, the equations to be solved are
(2.3)
+ =
(2.4)
+ =
(2.5)
= =
in which c is the speed of sound and K is the fluid bulk modulus (or the effective bulk
modulus in an elastic pipe). In this case, Eqs 2.3 & 2.4 can be developed into a typical
MOC form, giving two equations
(2.6a,b)
± = ±
(2.7a,b)
=±
respectively. This form of the MOC equations relates pressure changes to velocity
changes in much the same way that steady-state friction expressions can relate pressure
gradient to velocity. Note that the common step of replacing the pressure by pressure
head has not been taken here because it involves an ambiguity that is best left to a later
stage of the discussion.
For future reference, it is useful to declare that the justification for neglecting the
convective terms in the continuity equations is that they are small in typical water-
hammer applications. However, the reason for neglecting them is that they are
inconvenient in analyses using fixed-grid discretisation of space-time. We shall return to
this matter in Section 5. At this stage it is sufficient to record that, if the convective
terms are not neglected, Eqs 2.6a,b remain exactly as before, but the characteristic
directions become
(2.8a,b)
= ±
in which the suffices LA and RA denote averaging along the respective characteristics.
The particular form of averaging used is influenced by personal preference as well as by
historical publications and need not concern us here. It is the existence of the
inconvenient need for averaging that matters. This, together with uncertainty in the
values of parameters such as material properties and with approximations used to
represent phenomena contributing to the RHS of Eqs 2.1 and 2.2, is commonly (and
sensibly) used to justify simplifications that make our lives easier. Such simplifications
often include neglecting the convective terms and adjusting the assumed sound speed to
enable interpolation to be avoided.
ΔtLA
ΔtRA
L
R
Fig 2.1 Characteristic lines in x-t space
The advantages of avoiding interpolation are easily inferred from Fig 2.2, which
illustrates complications that its existence introduces. The most obvious is the additional
complexity of needing to estimate values of all relevant variables and parameters at the
points L and R when these do not coincide with grid points. In Fig 2.2a, the points Lx,Rx
imply space-line interpolation, the points Lt,Rt imply time-line interpolation, and the
points L*,R* denote grid points adjacent (spatially) to the solution point A. Fig 2.2b
shows loosely similar interpolation points for a diamond grid.
A A
L* R* L* R*
Lx Rx Lx Rx
Lt Lxt Rxt
Rt
In practical simulations, the MOC equations are usually integrated numerically over
finite intervals Δt. To do this, it is necessary to approximate the coefficient ρc, which,
strictly, is a variable. Often, as in acoustics, this coefficient is treated as a universal
constant, (ρc)0, say. If desired, however, it can be treated more accurately as quasi-
constant, based on local values of ρ and c at each location (x,t). Either way, we expose
ourselves to the layman’s curiosity as expressed above regarding the paradoxical
assumption of constant density in an analysis of waves in pipes.
The curiosity increases when authors take the further step of replacing the pressure p in
the compatibility equations by the pressure head h = p/(ρg) and express the outcome as
ℎ (3.1a,b)
± = ±
Strictly, h is a measure of the ratio p/ρ, not of p alone. Therefore, this substitution is not
strictly valid because it neglects an additional dρ/dt term. However, if we were obliged
to proceed in this rigorous manner, there would be no advantage in making the
substitution. Therefore we do not do so. Instead, we accept the ambiguity and enjoy the
practical benefits that the assumption of constant-density affords. That is, engineering
pragmatism determines our action. Even so, it is instructive to ask “How should density
be interpreted in the evaluation of the head?” Should we use local values of ρ or a fixed
value ρ0, say, independent of the pressure?
1 (3.2a,b)
± = ±
again applicable in the directions given by Eqs 2.7 or 2.8, as appropriate. These
equations are exactly equivalent to Eqs 2.6. Indeed, each can be derived directly from
the other using Eq 2.5. One inconvenience of the revised formulation is that it requires
our enquiring layman/student/whoever to regard density as constant even in an equation
in which it is one of the principal variables. Nevertheless, we know that this can be
justified so let us press on.
(3.3a,b)
= −2 ; = −2
1 (3.4a,b)
−2 ± = ±
This formulation has the big advantage that the coefficients of the differential terms are
simple constants and so the left hand sides of the MOC equations can be integrated
exactly over finite time intervals. In the special case of CRHS = 0 and MRHS = 0 (e.g.
simple inviscid flow), no approximation whatsoever is required to integrate the
equations.
One technical reason is that engineers are usually more interested in determining
variations in pressure than variations in the sound speed. As a consequence, the sound-
speed formulation requires an additional step at the output stage, namely converting from
c to p. This step is not difficult, but it is an irritation, especially because the p = p{c}
relationship is logarithmic, namely dp = 2Kd(ln(c)). This can be integrated exactly
between a fixed reference state and the current state. Notice, however, that this is only an
output matter. The derived value of pressure is not used in the ensuing calculation
process at other (x,t) so it does not influence the solution accuracy.
A third technical reason is that the (c,U) formulation can be rather inconvenient at
boundaries, especially those that are prescribed in terms of pressure. At such boundaries,
the above relationship between pressure and sound speed must be used in the solution
process itself. Furthermore, this can be done exactly only in an iterative manner. Also,
even if the p = p{c} relationship is linearized to make an explicit solution possible, as
would often be reasonable, the need for the additional equation is nevertheless an
unwelcome complication. Moreover, further complications arise in boundary conditions
that involve the density. If this is treated as a variable (necessary if an “exact” solution is
sought), it must be substituted by K/c2, thereby usually removing the possibility of a truly
linear boundary condition.
Our final technical reason is by far the most important. We have deliberately left it to the
end in the hope of reminding any unwary reader how easy it is to overlook the obvious
(as viewed retrospectively). The second line of this Section (3.2) includes the words “if
the bulk modulus is constant”. This turns out to be a rather big “if” that is exposed
starkly in Eqs 3.3 which imply that the sound speed decreases with increasing pressure
and with increasing density. This is contrary to the behaviour of common liquids,
including water. That is, it is physically unrealistic to treat the bulk modulus as constant
in the development of dp/dc and dρ/dc even though this approximation is usually
acceptable in the formulation of dp/dρ given in Eq 2.5. We leave the reader to figure out
why is this so. The reason is obvious with hindsight, but the process of discovering it
will bring valuable benefit in its own right to anyone for whom the reason is not already
clear [Hint: study the implications of Eq 2.5].
Before moving to the next topic, it is worth noting that all of the above issues (except
pressure head) also arise in “compressible” fluid-hammer analyses. Typically, the
equations for such flows are more complex than those describing “incompressible”
water-hammer. However, there is a notable exception, namely waves in homentropic
flows of a perfect gas (N.B.: “homentropic” implies not only that changes occur
isentropically, but also that all fluid in the system has the same specific entropy). For
this special case, the (c,U) formulation of the equations is almost the same as Eqs 3.4, the
only difference being the replacement of the coefficient “–2” by a positive constant. In
the particular case of a gas for which the ratio of specific heats may be taken as 1.4 (a
common assumption for air), the value of the constant is “+5”. This form of the MOC
equations is quite convenient in gas flow simulations because the alternative formulations
in (p,U) or (ρ,U) lose the simplifying benefits that they have in “incompressible” water-
hammer.
The focus of Section 3 is on LHS terms in the compatibility equations. Let us now look
at the RHS terms, namely the non-differential ones. These terms exist to describe
physical phenomena such as friction, source/sink flows and spatial area change. When
they are small in comparison with the LHS terms, their inclusion in the numerical
integration is straightforward. When they are large, however, they can be troublesome,
although it is often possible to overcome such problems by, for example, choosing
numerical approximations thoughtfully and/or limiting the integration time steps.
Herein, we are not concerned primarily with matters of stability, etc. Instead, we have a
new little stone to chuck into the pond. Its ripples might provide amusement for some.
Many RHS terms can be expressed in the generic form (using the forward characteristic
as an exemplar):
+ + | | (4.1)
where the coefficients αi can be reduced to numerical constants. This is commonly the
case in practical computations although it is often achieved only with the help of some
approximations. Possible physical interpretations of the three terms include body forces,
laminar friction and turbulent friction, but the following discussion is not restricted to
these particular interpretations. Furthermore, in everything that follows, the
methodology applied to |U|U terms could equally well be used for U2 terms if these also
exist. Indeed, a similar methodology can also be used for higher-order terms (e.g.
U3 = U2.U).
+ + | | (4.2)
+ + | | (4.3)
+½ ( + )+½ (| | +| | ) (4.4)
The first two of these are end differences, one using values that are already known and
one using unknown values at the current solution point. They may be expected to be of
relatively low order accuracy, but that objection can often be overcome by using
sufficiently small integration steps. When developing code, the second (Eq 4.3) requires
a little more effort than the first, especially if the |U|U term exists - because the solution
+ + | | (4.5)
+½ ( + )+ | | (4.6)
+ (4.7)
This method of linearising higher order terms – and hence avoiding the need for
iterations – provides stability and also has some of the character of a geometrical central
difference approximation. It is the method on which the remainder of this section
focusses. Note, however, that similar general outcomes can be obtained with other
possible methods of approximating RHS terms.
+ [( ) − ] = + [( ) + ] (4.10)
− [( ) + ] = − [( ) − ] (4.11)
Numerically, Eqs 4.10 and 4.11 are treated as linear simultaneous equations in pA and UA.
If the products ρc and the coefficients βi are treated as constants, the solution is direct.
Alternatively, if ρc or βi are treated as variable, the process will be iterative, although it
may usually be expected to converge rapidly.
( − ) + [( ) − ]( − )=[ + ] (4.13)
Now compare this with the equation that would have been obtained if the RHS terms had
been developed using end difference methods typified by 4.2 instead of 4.5, namely
That is, although additional terms involving UA–UL were developed as RHS terms in the
numerical equations, they will behave as if they had been developed from LHS terms in
the analytical equations. In other words, if ρc and the coefficients βi are regarded as
constants, equations identical in all respects to 4.10 and 4.11 would have been obtained if
the original compatibility equation 2.6a had been of the form:
(4.15)
+ [( ) − ] =[ + ]
For completeness, we note in passing that, in the limit as Δt tends to zero, this equation
tends to the original differential equation. The existence of the numerical product βΔt on
the LHS looks strange, but it is a simple constant so it is easily handled mathematically.
Nevertheless, up to this point, its existence arises solely because of the numerical
integration, not because of physics.
We now move to the crux of the matter by asking what change would need to be made to
the original equations 2.1 & 2.2 (i.e. to the representation of the physics) to cause Eq 4.15
to arise. Clearly, the product ρc would need to change and this implies a change in the
bulk modulus and hence the sound speed c (NB: A change in the density alone would not
have the desired effect because it would also influence c through Eq 2.5). Specifically,
we would replace the true sound speed by c – β1LΔt/ρ. That is, the possible re-writing of
RHS terms such that they depend upon changes in velocity has consequences that change
the MOC compatibility equations in a manner that has the same effect as a physical
change in sound speed – and hence in wave speed.
The authors are not suggesting that the implied change in sound speed is physically real.
On the contrary, it clearly is not. Furthermore, the overall solutions from the two
analyses over multiple time steps would be different – because the gradients of the
characteristics would change if the modified wave speed were used whereas they are not
changed in the above development. Nevertheless, it is widely recognised that changes in
the MOC compatibility equations have greater influence on overall solutions than
changes in the gradients of the characteristics. Therefore, it is reasonable to consider the
existence of “UA” in the RHS terms as almost equivalent to a change in sound speed. For
completeness, it is noted that the existence of “pA” in the RHS terms would have a
similar effect because pA can be interpreted as pL + (pA – pL).
Second, we declare a further reason why the effect should not always be interpreted in
general as a simple change in the effective sound speed. This is because the two
compatibility equations (LA & RA) imply different adjustments, namely c – β1LΔt/ρ and
c + β1RΔtRA/ρ. The two adjustments will be equivalent only when β1R = β1L. This will
be true for some physical causes of RHS terms, but not for others. Typically, it may be
approximately true for terms such as friction that “belong” to the momentum equation,
but not for terms such as mass sources/sinks that also appear in the continuity equation.
Attention so far has focussed on the MOC compatibility equations. We now ask whether
the equations for the characteristic directions have any surprises in store. To pursue this
question, we examine the implications of the simplification made to approximate the
“true” characteristic directions dx/dt = U±c in Eq 2.8 by the constant directions
dx/dt = ±c in Eq 2.7.
Figure 5.1 depicts characteristics passing through grid points L,R in a rectangular grid
chosen so that Δx/Δt = c. Both characteristics will pass through the intermediate grid
point Δt later if, and only if, their average gradients satisfy (dx/dt)LA = (dx/dt)RA = c and
this will rarely be true unless their gradients are everywhere equal to c, as implied by
Eq 2.7. In all other cases, the characteristics will meet at some other point, as illustrated
in Fig 5.1b. Furthermore, the characteristics will usually be curved, but this will be of
secondary importance in practical simulations of low Mach number flows so we shall
focus herein only on the average gradients.
L R L R
In practical simulations that include the more general form of the characteristic
directions, it is usual to make an adjustment so that the characteristics do pass through
the point A. This can be achieved by reducing the integration time step so that Δt ≤ Δx/c
everywhere, typically yielding one of the geometries depicted in Fig 5.2, which illustrates
space-line and time-line interpolation in a rectangular grid. Similar outcomes obtain
when diamond grids are used. The need for interpolation is unavoidable when
dx/dt = U±c is used instead of dx/dt = c, and this can be a significant computational
penalty in cases where it is important for the computer software to be computationally
fast. Also, it complicates the coding itself, albeit not necessary greatly. However, there
are other possible causes of a need for interpolation – e.g. (i) different values of Δx in
different pipes in a network or (ii) different effective sound speeds when allowance is
made for different pipe elasticities or restraints. When such effects exist, the penalty for
using dx/dt = U±c instead of dx/dt = c is small.
A A
Lx Rx
Lt
Rt
(a) Space-line interpolation (b) Time-line interpolation
Fig 5.2 Use of interpolation to ensure intersection at a grid point A
In simple cases such as simulations for single pipes, interpolation can be avoided
altogether by using irregularly spaced solution points such as A as the starting points for
subsequent steps, thereby allowing a natural grid to develop. However, this is not
possible in general pipe networks because the independent natural grids in adjacent pipes
at any particular internal boundary do not conveniently intersect the boundary at the same
instant. Moreover, even in the cases of a single pipe, output will usually be required at
(5.1a,b)
± =0
irrespective of the formulation used for the characteristic directions. However, Fig 5.1
shows that the solution of these two equations will be deemed to apply at different points
in space-time, depending upon the choice of dx/dt. It follows that an inevitable
consequence of using the simplified gradient is to displace the solution in both space and
time.
Reverting to the more general case when the RHS of one or both of the continuity and
momentum equations is non-zero, the outcome is a little more complicated because the
different values of, say, (tA–tL) seen in Fig 5.1 a & b will cause differences in the
calculated values of pA,UA as well as differences in the locations of the solution points
themselves. That is the use of the approximate characteristic will not only displace the
solution, but will also distort its amplitude.
Now consider the more practical numerical grid usage depicted in Fig 5.2 or in
equivalent diamond grid applications. If space-line interpolation is used, the integration
time steps will be predetermined and so there will no longer be a difference in the time
steps used in the two methods. However, the elimination of this difference is achieved
only at the expense of accepting non-grid locations of the point L and so causing the
integrations to be undertaken with interpolated values of pL,UL (and likewise for the other
characteristic). Hence the numerical integrations of the RHS terms will again be
different and so the calculated values of pA,UA will again be different.
Water-hammer specialists might justifiably argue that the effect we are discussing is
likely to be too small to be of concern in typical practical applications and, indeed, even
in some more academic ones. We would agree with this argument. However, we would
also caution against allowing the conclusion to influence more general applications of
MOC methods. For example, the effect can be important in quite low speed gas flows.
The first author has a special interest in unsteady airflows in railway tunnels. One of the
most important phenomena encountered with high speed trains in long tunnels would be
almost totally suppressed in an analysis based on dx/dt = ±c even though the associated
Mach number is usually smaller than 0.02.
5.2 Steady-flows
Many discussions over the years have drawn attention to a perceived inability of MOC
solutions to reproduce steady state conditions exactly. It is widely recognised that, if
everything were modelled correctly, true steady-state should be the asymptotic result of
independent solutions with successively smaller grid sizes. However, this outcome is
commonly not obtained when the simulations are undertaken. Let us explore some
possible reasons for this.
Possibility-2: The method of integrating the RHS terms in the MOC equations is invalid.
This could be the case if the method is chosen without due care, but any reasonable
method should be acceptable when sufficiently small grid sizes are used.
Possibilities 1-3 collectively include all possible sources of MOC error that the authors
can identify and yet they do not provide an explanation for differences between the MOC
steady state solution and other methods of solution (unless the stated precautions are not
taken). As a consequence, we must look elsewhere for the error – and there is only one
place to look, namely the solution with which the MOC outcome is being compared.
Of course, this conclusion is already well known to specialists and so is the most likely
source of error in a typical (non-MOC) steady state solution. Quite simply, it is an
expectation that the velocity will be uniform in a steady flow in a uniform pipe. In
reality, the pressure will vary along the pipe (because of friction and elevation change,
say). As a consequence, the density will also vary and, to satisfy continuity, the velocity
will vary inversely with the density. In liquid flows, the velocity changes will be small,
but they will not be zero. That is, ‘steady-state’ flow in a uniform pipe is not equivalent
to ‘uniform velocity’.
6 CLOSING REMARKS
Properties of the method of characteristics have been assessed from a somewhat unusual
standpoint. In part, the purpose has been to investigate why the method is so successful
in the simulation of liquid flows in pipelines. In part, it has been to investigate what it
can and cannot do. The key outcomes may be summarised as follows.
• The formulation of the equations with pressure and velocity (p,U) as the principal
variables is more accurate than alternative formulations that, at first sight, appear to
be equivalent. The alternative formulations that have been considered are (h,U),
(ρ,U) and (c,U).
• When the RHS terms are evaluated in a manner that includes the velocity at the
solution point, the outcome is equivalent to a change of wave speed in the
compatibility equations. Strictly, this should be allowed for in the evaluation of
characteristic directions – although the authors have not done this themselves and
are not aware of any instances where others have done it.
ACKNOWLEDGEMENT
The authors’ desire to explore MOC in the manner presented above is a consequence of
discussions over a long timescale (decades) with students, colleagues, fellow researchers
and conference delegates. We give special thanks to all who stubbornly refused to accept
anything we said at face value.
REFERENCES
The authors have taken an active decision not to cite any references in the body of the
paper. This is because we seek to provoke discussion, not to provide a comprehensive
(or even a selective) review of previous work. However, we should draw attention to the
following discussions on the ability of MOC to model steady flow accurately.