Jtam.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

49, 1, pp. 135-158, Warsaw 2011

IMPROVED METHOD FOR SIMULATING TRANSIENTS OF

TURBULENT PIPE FLOW

Zbigniew Zarzycki

Sylwester Kudźma

Kamil Urbanowicz

West Pomeranian University of Technology, Faculty of Mechanical Engineering and Mechatronics,

Szczecin, Poland; e-mail: zbigniew.zarzycki@zut.edu.pl; sylwester.kudzma@zut.edu.pl;

kamil.urbanowicz@zut.edu.pl

The paper presents the problem of modelling and simulation of tran-
sients during turbulent fluid flow in hydraulic pipes. The instantaneous
wall shear stress on a pipe wall is presented in the form of integral co-
nvolution of a weighting function and local acceleration of the liquid.
This weighting function depends on the dimensionless time and Rey-
nolds number. Its original, very complicated mathematical structure is
approximated to a simpler formwhich is useful for practical engineering
calculations.The paper presents an efficientway to solve the integral co-
nvolution based on themethod given by Trikha (1975) for laminar flow.
An application of an improved method with the use of the Method of
Characteristic for the caseof unsteadyflow(water hammer) is presented.
This method is characterised by high efficiency compared to traditional
numerical schemes.

Key words: unsteady pipe flow, transients, waterhammer, efficient nu-
merical simulation

Notation

c0 [ms
−1] – acoustic wave speed

p [kgm−1s−2] – pressure

s [s−1] – Laplace operator

t [s] – time

v [ms−1] – instantaneous mean flow velocity

t̂= νt/R2 [–] – dimensionless time



136 Z. Zarzycki et al.

v [ms−1] – average value of velocity in cross-section of pipe

vz [ms
−1] – axial velocity

w [–] – weighting function

z [m] – distance along pipe axis

L [m] – pipe length

R [m] – radius of pipe

Re=2Rv/ν [–] – Reynolds number

µ [kgm−1s−1] – dynamic viscosity

ν [m2s−1] – kinematic viscosity

λ [–] – Darcy-Weisbach friction coefficient

ρ0 [kgm
−3] – fluid density (constant)

τw,τwq,τwu [kgm
−1s−2] – wall shear stress, wall shear stress for quasi-steady

flow, unsteady wall shear stress, respectively

Ω=ωR2/ν [–] – dimensionless frequency

1. Introduction

There are many problems concerning the issue of unsteady flow in pipes, and
they include such phenomena as: energy dissipation, fluid-structure interac-
tions, cavitation (for example column separation) and many others. It would
bevery difficult, if not impossible, to account for all the phenomena in just one
mathematical model. This paper is devoted to the problem of energy dissipa-
tion and it concerns unsteady friction modelling of the liquid in pipes during
turbulent liquid pipe flow.

The commonly used quasi-steady, one-dimensional (1D) model in which
the wall shear stress is approximated bymeans of Darcy-Weisbach formula is
correct for slow changes of the liquid velocity field in the pipe cross-section, i.e.
either for small accelerations or for low frequencies (if the flow is pulsating).
The formula is also correct for the waterhammer effect for its first wave cycle.
Thismodel does not take into account the gradient of speed changing in time,
which is directed into the radial direction, and so, as a result, the changing
energy dissipation is not taken into account either.

Models which represent unsteady friction losses can be divided into two
groups. The first are those inwhich the shear stress connectedwith the unste-
ady friction term is proportional to the instantaneous local acceleration (Daily
et al., 1956; Carrtens and Roller, 1959; Sawfat and Polder, 1973). This group



Improved method for simulating transients... 137

contains a popular model developed by Brunone et al. (1994, 1995) in which
the unsteady friction term is a sum of two terms: the first of which is pro-
portional to the instantaneous local acceleration and the other is proportional
to the instantaneous convective acceleration. Vitkovsky et al. (2000) improved
Brunone’s model by introducing a sign to the convective acceleration.

In the second group of models, the unsteady wall shear stress depends on
the past flowaccelerations. Thesemodels are based on a two-dimensional (2D)
equation of motion and they take into account the velocity field which can be
changeable in time in the cross-section of the pipe. Zielke (1968) proposed
for the first time a laminar flowmodel in which the instantaneous shear stress
dependson the instantaneouspipeflowacceleration andtheweighting function
of the past velocity changes. In the case of unsteady turbulent flow,models of
friction losses aremainly based on the eddy viscosity distribution in the cross-
section of the pipe determined for a steady flow. The followingmodels should
be listed: two and three region models (Vardy, 1980; Vardy and Brown, 1995,
1996, 2003, 2004, 2007; Vardy et al. 1993) and as well as Zarzycki’s four-layer
model (Zarzycki, 1993, 1994, 2000). Lately, Vardy and Brown (2004) used an
idealised viscosity distribution model for a rough pipe flow.

Unsteady friction models are used for simulation of transients in pipes.
The traditional numerical method of unsteady friction calculation proposed
by Zielke (1968), requires a large amount of computer memory and is time
consuming.Therefore,Trikha(1975) andthenSuzuki (1991) andSchohl (1993)
improved Zielke’s methodmaking it more effective.

In the case of unsteady turbulent flow, Rahl and Barlamond (1996) and
thenGhidaoui andMansour (2002) usedVardy’smodel to simulate thewater-
hammer effect. These methods provide only approximated results. They are
not very effective and cannot be considered to be satisfactory.

Theaimof thepresentpaper is todevelopaneffectivemethodof simulating
shear stress for unsteady turbulent pipe flow, amethod that would be similar
to that proposed by Trikha and Schoohl for laminar flow. That also involves
developing an approximated model of the weighting function, which can be
useful for effective simulation of transients.

2. Equations of unsteady turbulent flow

Unsteady, axisymmetric turbulent flow of a Newtonian fluid in a long pipe
with a constant internal radius and rigid walls is considered. In addition, the
following assumptions are made:



138 Z. Zarzycki et al.

• constant distribution of pressure in the pipe cross-section,
• body forces and thermal effects are negligible,
• mean velocity in the pipe cross-section is considerably smaller than the
speed of sound in the fluid.

The following approximate equation, presented in the previous papers by
Ohmi et al. (1980–1981, 1985), Zarzycki (1993), which omits small terms in
the fundamental equation for unsteady flow and in which the Reynolds stress
term is described by using eddy viscosity νt, is used

∂vz
∂t
=− 1
ρ0

∂p

∂z
+νΣ
∂2vz
∂r2
+
(
r
∂νΣ
∂r
+
νΣ
r

)∂vz
∂r

(2.1)

where vz, p – averaged in time, respectively: velocity component in the axial
direction and pressure, ρ0 – density of the liquid (constant), t – time, z –
distance along the pipe axis, r – radial distance from the pipe axis.

The modified eddy viscosity coefficient νΣ is the sum of the kinematic
viscosity coefficient ν and the eddy viscosity coefficient νt, that is

νΣ = ν+νt (2.2)

It is assumed that changes of the coefficient νt resulting from the flow change
in time are negligible, which means that νt is only a function of the radial
distance r. In practice, this limits the analysis to relatively short periods of
time.

2.1. Turbulent kinematics viscosity distributions

In order to describe the distributionof the turbulentviscosity coefficient νt
throughout the pipe cross-section, the flow region is divided into several layers
and, for each of them, the function νt(r) has to be determined. By analogy
to steady pipe flow, the following regions of flow can be distinguished: VS
(Viscous Sublayer), BL (Buffer Layer), DTL (Developed Turbulent Layer)
and TC (Turbulent Core) (Hussain and Reynolds (1975)). This is shown in
Fig.1. The distribution of eddy viscosity presented in Fig.1 was created on
the basis of experimental data of the coefficient νt provided by Laufer (for
Re=4.1·105) andNunner (for Re=3·104) and publishedbyHinze (1953) as
well as data about a region close to thewall, whichwas provided bySchubauer
and published by Hussain and Reynolds (1975). That was presented in depth
by Zarzycki (1993, 1994).



Improved method for simulating transients... 139

Fig. 1. Schematic representation of the four-regionmodel of νt

The radial distances r1, r2 and r3 from the pipe centre line are as follows

rj =R−y(y+j ) j=1,2,3
y+1 =0.2R

+ y+2 =35 y
+
3 =5

(2.3)

where: y+ = yv∗/ν is the dimensionless distance y, y – distance from the
pipewall, v∗ =

√
τws/ρ0 – dynamic velocity, τws –wall shear stress for steady

flow, R+ =Rv∗/ν – dimensionless radius R of the pipe.

The coefficient of turbulent viscosity νt is expressed for particular layers
(except the buffer layer) as follows

νt =αir+βi i=1,3,4

α1 =0 α3 =−
νt1−νt2
r2−r1

α4 =0

β1 =0 β2 = νt2−α3r2 β4 = νt1

(2.4)

For the buffer layer (r2 ¬ r < r3), νt is expressed as follows

νt =0.01(y
+)2ν (2.5)

Expressions, which define the quantities νt1, νt2, r1, r2 and r3 have the follo-
wing forms

νt1 =0.016
v̂∗√
2
ν νt2 =12.25ν (2.6)

and

r1 =0.8 r2 =
(
1−
140
√
2

v̂∗

)
R r3 =

(
1−
20
√
2

v̂∗

)
R (2.7)



140 Z. Zarzycki et al.

Thequantities νt1,r2 and r3,whichdefinethedistributionof the coefficient νt,
depend on a dimensionless dynamic velocity v̂∗, which is determined in the
following way

v̂∗ =
4
√
2R

ν
v∗ =

√
λsRes (2.8)

where: λs is the friction coefficient for steady flow, Res =2Rvs/ν – Reynolds
number for steady flow, vs – mean fluid velocity in the pipe cross-section.

Since we have λs = λs(Res) then, for smooth pipes, the profile of the
coefficient νt in the pipe cross-section depends only on the Res number.

However, if – instead of the Res – we assume the instantaneous Reynolds
number Re = 2Rv/ν (v – instantaneous, averaged in the pipe cross-section
fluid velocity) and – instead of λs – we assume a quasi-steady friction coef-
ficient λq = λq(Re), then the presented distribution of the coefficient νt will
become a quasi-steady one. In this paper λq is determined from the Prandtl
formula

1
√
λq
=0.869ln(Re

√
λq)−0.8 (2.9)

2.2. Momentum equations for particular layers

For the four-layer model of the flow region, equation (2.1) breaks up into
four equations. Introducing the differential operator D = ∂/∂t into equation
(2.1), we obtain:

— for the viscous sublayer

Dvz1 =−
1

ρ0

∂p

∂z
+ν
(∂2vz1
∂r2

+
1

r

∂vz1
∂r

)
(2.10)

— for the remaining layers of the flow region

Dvzi =−
1

ρ0

∂p

∂z
+νΣ
∂2vzi
∂r2
+
(∂νΣ
∂r
+
νΣ
r

)∂vzi
∂r

(2.11)

where: i = 2,3,4 (i = 2 for the BL, i = 3 for the DTL and i = 4 for the
TC), vzi is the axial component of the velocity for particular layers of the flow
region, νΣ = ν+νt – effective coefficient of turbulent viscosity.

Equations (2.10) and (2.11) (for i=3,4) canbe resolvedusing relationship
(2.4) into the modified Bessel equations. Equation (2.11) (for i = 2), taking
into account the fact that



Improved method for simulating transients... 141

r

R−r
≫ 1
2

for r2 ¬ r < r3 (2.12)

assumes the form of a differential equation of Euler type.

The solutions to the momentum equations for particular layers are shown
in papers (Zarzycki, 1993, 1994, 2000).

3. Wall shear stress and weighting function

Integrating equation (2.1) over the pipe cross-section (from r =0 to r =R)
and considering that ∂/∂t=D, we obtain

ρ0Dv+
∂p

∂z
+
2

R
τw =0 (3.1)

By determining the wall shear stress at the pipe wall τw as follows

τw =−ρ0ν
∂vz1
∂r

∣∣∣∣
r=R

(3.2)

we obtain a relationship of the following type

τw = f(D,Re)
∂p

∂z
(3.3)

Equations (3.1) and (3.3) enable one to determine transfer function Ĝτv that
describes hydraulic resistance, relating the pipewall shear stress and themean
velocity

Ĝτv(D̂,Re)=
τ̂w
v̂

(3.4)

where v̂ = (R/ν)v, τ̂w = [R
2/(ρ0ν

2)]τw, p̂= [R
2/(ρ0ν

2)]p denote dimension-
less velocity, wall shear stress, and pressure, respectively.

A precise mathematical form of the transfer function Ĝτv is very complex
and was presented by Zarzycki (1994).

The quantity D̂=(R2/ν)∂/∂t presents the dimensionless differential ope-
rator.

Taking into account transfer function (3.4) for D̂ = ŝ (where ŝ= sR2/ν
is the dimensionless Laplace operator) one can express the stress τw as a
sum of quasi-steadywall shear stress τwq and integral convolution of the fluid
acceleration ∂v/∂t and the weighting function w(t), i.e. as a convolution of



142 Z. Zarzycki et al.

thefluid acceleration ∂v/∂t and the step response g(t) (Zielke, 1968; Zarzycki,
1993, 2000), i.e.

τw(t)= τwq+
2µ

R

t∫

0

∂v

∂t
(u)w(t−u) du (3.5)

where τwq is described by the following expression

τwq =
1

8
λqρ0v|v| (3.6)

The graph of weighting function w(t) for different Reynolds numbers and for
the laminar flow is presented in Fig.2.

Fig. 2. The weighting functions for turbulent and laminar flow

The weighting function approaches zero when time tends to infinity. This
tendency becomes stronger for higher Reynolds numbers.

4. Approximated weighting function model

The weighting function has a complex mathematical form and it is not ve-
ry useful to calculate and simulate unsteady turbulent flow (Zarzycki, 1994,
2000). Zarzycki (1994) approximated it to the following form



Improved method for simulating transients... 143

wapr =C
1√
t̂
Ren (4.1)

whereC =0.299635, n=−0.005535.
However, numerical calculations that include this function lead to some

complications (dividing by zero) as time and the Reynolds number are to
be found in the denominator. Hence, in order to improve the effectiveness of
calculations, the function will be approximated with a simpler form. Additio-
nally, it is suggested that the approximated function shouldmeet the following
assumptions:

1) its form should be similar toTrikha’s (Trikha, 1975) or Schohl’s (Schohl,
1993)models (methods of simulating shear stress based on thesemodels
are very effective),

2) its mathematical form should allow an easy way of carrying out Laplace
transformation. This is important for simulating oscillating or pulsating
flow by thematrix method (Wylie and Streeter, 1993).

In order to meet the above presented criteria, a number of possible forms
of the functionwere investigated. The formwhichwas eventually selected can
be given by

wN(t̂)= (c1Re
|c2|+ c3)

n∑

i=1

Aie
−|bi|̂t (4.2)

where i=1,2, . . . ,n.

Determining the coefficients of the weighting function was carried out in
two stages in the following way:

a) Thefirst stage. In the rangeof theReynolds numbersRe∈ 〈2·103,107〉,
a several equally distributed points were chosen, and for each of them
weighting functionswere estimated usingStatistica programme.Estima-
tion was conducted for the range of dimensional time t̂∈ 〈10−6,10〉.

b) The second stage. Two-dimensional weighting function w(t̂,Re) was
estimated usingMatlab’s module for optimization. For the optimization
procedure, the following assumption was formulated:

– to determine the decision variables vector x = [x1,x2, . . . ,x19],
where: x1 = c1, x2 = c2, x3 = c3, x4 = A1, x5 = b1, x6 = A2,
x7 = b2, . . . , x18 =A8, x19 = b8 (are coefficients of (4.2))

– which minimizes the vector objective function

minf(x)= [f(x)] (4.3)



144 Z. Zarzycki et al.

where

f(x)=
r∑

i=1

(wN(i)(t̂,Re)−wi(t̂,Re)
wi(t̂,Re)

)2
(4.4)

and satisfies the following limitation

g(x) =

10
7∫

2·103

[ 1∫

10−5

wN(t̂,Re)dt̂

]
dRe−

10
7∫

2·103

[ 1∫

10−5

w(t̂,Re)dt̂

]
dRe=0

(4.5)

In relation (4.4) the arguments t̂ andRewere selected so that the determined
points (pairs ofnumbers)on theplane t̂−Rewoulduniformlycover the interval
t̂∈ 〈10−5,10−1〉 and Re= 〈2·103,107〉, and that theywould be different from
the points selected at the first stage. A total of r=1000 points were used in
the calculations. The values of coefficients in Eq. (4.2), determined at the first
stage, were used as the starting points (the values of elements of the decision
variables vector) for the optimization.

Once the optimization problemswere formulated in thewaydescribed abo-
ve, calculations were conducted usingMatlab and its function ”fgoalattain”.

Details of the whole process of a new weighting function estimation were
described by Zarzycki and Kudźma (2004).

On the basis of the above assumptions and calculations that were conduc-
ted later, the final form of the weighting function was established

wN(t̂)= (c1Re
c2 + c3)

8∑

i=1

Aie
−bit̂ (4.6)

c1 = −13.27813; c2 = 0.000391; c3 = 14.27658; A1 = 0.224, A2 = 1.644,
A3 =2.934, A4 =5.794, A5 =11.28, A6 =19.909, A7 =34.869, A8 =63.668;
b1 = 0.10634, b2 = 8.44, b3 = 88.02, b4 = 480.5, b5 = 2162, b6 = 8425,
b7 =29250, b8 =96940.

Thepresent formof theweighting function contrary to the oneput forward
in (Zarzycki and Kudźma, 2004) contains eight terms. As it turned out, the
effective method of calculating wall shear stress presented later in this work is
more sensitive to inaccurate fitbetween theproposed functionand theoriginal.
Therefore, in order to ensure correct mapping, the number of elements of the
new weighting function was increased from six to eight.

In addition, a comparison was made of new expression (4.6) course and
Zarzycki’s model with the weighting function of Vardy and Brown (2003)



Improved method for simulating transients... 145

w(t̂)=
1

2
√
πt̂
exp
(
−
t̂

C∗

)
(4.7)

where:C∗ =12.86/Reκ; κ= log10(15.29/Re
0.0567).

The accuracy of Zarzycki’s mapping of the weighting function by the new
function is as follow:

• For a range of Reynolds numbers between 2 · 103 ¬ Re ¬ 107 and
dimensionless time between 10−5 ¬ t̂¬ 10−1, the relative error satisfies
condition

∣∣∣
wN(t̂)−w(t̂)
w(t̂)

∣∣∣¬ 0.05

Graphical comparison of the newly estimated function andmodels of weighing
function by Zarzycki and Vardy-Brown are presented in Fig.3. The relative
error of the new function and Zarzycki’s model is presented in Fig.4. As it is
shown, the error in some areas is up to 5%but it has to bementioned that the
run of the new function is between the models of Vardy-Brown and Zarzycki
weighting functions. There is no clear evidence of better consistency with the
experiment for awide range ofReynoldsnumbersup tonow, so it is ambiguous
which models should be mapped more precisely by the new function. The
authors did not want to complicate the form of weighting expression through
an increase of the number of exponential components thus left it in form (4.6).
But, it is worth tomention that it is possible to increase accuracy ofmapping
the original weight in a wider range of Reynolds numbers or dimensionless
time in the case of newmore precise and experimentally tested models.

The lowest limit of usage of the new weighting function for dimensionless
time is t̂ = 10−5. For the presented method of transients simulations it is a
satisfactory value from the point of engineering view. It permit simulations
of low viscosity fluids in quite large pipe radius (e.g. for fluid which viscosity
ν =1mm2/s acceptable diameter = 0.3m). It has to be also mentioned that
such a big radius of the pipeline has usually long length, what means, in
consequence, low frequency of transients pressure oscillations. It allows one to
take a quite big time step in numerical calculations, which is in the acceptable
range of dimensionless time.

In order to quantitatively assess a mapping of Zarzycki’s model by the
new expression for selected Reynolds numbers, a relative error is calculated
according to the expression

relative error=
∣∣∣
wN(t̂)−w(t̂)
w(t̂)

∣∣∣ (4.8)



146 Z. Zarzycki et al.

Fig. 3. Comparison of Vardy-Brown’s, Zarzycki’s and new weighting function runs
for different Reynolds numbers: (a) Re=104, (b) Re=106

Fig. 4. The relative error between the new function and Zarzycki’s weighting
function model for Re=104

Results of its usage are shown for the case of Reynolds number Re = 104 in
Fig.4.

In the case of harmonically pulsating flow, a transfermatrixmethod is usu-
ally used (Wylie and Streeter, 1993). For this method, an adequate (laminar
or turbulent) weighting function transformed into Laplace domain is needed.



Improved method for simulating transients... 147

Next, the Laplace variable is substitutet by ŝ= jΩ.

Conducting the abovementioned transformations, the newweighting func-
tion takes the following form

wN(jΩ)= (c1Re
c2 + c3)

8∑

i=1

Ai
jΩ+ bi

(4.9)

where: Ω=ωR2/ν is the dimensionless angular frequency.

In Fig.5, a comparison of real and imaginary parts of newweighting func-
tion (4.9) with Zarzycki and Vardy-Brown models is presented.

Fig. 5. Comparison of the new weighting function with models by Zarzycki and
Vardy-Brown for Re=104

To sumup, it should bepointed out that newapproximating function (5.3)
well reflects the weight of the unsteady friction model for turbulent flow for
t̂ ­ 10−5 and 2000 ¬ Re ¬ 107. However, as far as frequency is concerned,
the presented relation is right for a dimensionless frequency Ω ­ 10. These
particular ranges are most important as far as engineering applications and
precision of numerical simulations are concerned.

5. Improvement of the approximation method for simulations of

wall shear stress

Numerical calculations of the unsteadywall shear stress (convolution of the lo-
cal liquid acceleration and theweighting function) in the traditional approach
is both time-consuming and takes a lot of computermemory space (Suzuki et
al., 1991). InpapersbyTrikha (1975) andSchohl (1993) another,more efficient
method of calculating unsteady shear stress for laminar flowwas presented.



148 Z. Zarzycki et al.

This section presents the implementation of Trikha’s method into calcula-
tions of wall shear stress for turbulent flow using the developed here form of
weighting function (4.6).
Weighting function (4.6) can be given by

wN(t)=
8∑

i=1

wi (5.1)

where wi =(c1Re
c2 + c3)Aie

−biνt/R
2

.
Using the above form of theweighting function, the unsteady term of pipe

wall shear stress τwn can be written as

τwu =
2µ

R

8∑

i=1

yi(t) (5.2)

where

yi(t)= (c1Re
c2 +c3)Aie

−bi
ν

R2
t

t∫

0

e
bi
ν

R2
u∂v

∂u
(u) du (5.3)

The Reynolds number was used in the above expressions as a variable to be
updated at every iterative step. This approach was also suggested by Vardy
et al. (1993).
Because in the characteristicsmethod calculations are carried out for some

fixed time intervals, then

yi(t+∆t)=

t+∆t∫

0

wi(t+∆t−u)
∂v

∂u
(u) du=

=e
−bi

ν

R2
∆t
(c1Re

c2 + c3)Aie
−bi

ν

R2
t

t∫

0

e
bi
ν

R2
u∂v

∂u
(u) du

︸ ︷︷ ︸
yi(t)

+ (5.4)

+(c1Re
c2 + c3)Aie

−bi
ν

R2
(t+∆t)

t+∆t∫

t

e
bi
ν

R2
u∂v

∂u
(u) du

︸ ︷︷ ︸
∆yi(t)

While comparing relations (5.4) and (5.3), we obtain

yi(t+∆t)= yi(t)e
−bi

ν

R2
∆t
+∆yi(t) (5.5)



Improved method for simulating transients... 149

where

∆yi(t)= (c1Re
c2 +c3)Aie

−bi
ν

R2
(t+∆t)

t+∆t∫

t

e
bi
ν

R2
u∂v

∂u
(u) du (5.6)

If in the above given expression describing ∆yi(t), we assume that v(u) is
a linear function [v(u) = au+ b] in the interval from t to t+∆t, then
its derivative calculated after time given by ∂v(u)/∂u can be treated as
a constant whose value can be determined from the following expression:
[v(t+∆t)−v(t)]/∆t. On the strength of this assumption, we can write

∆yi(t)≈ (c1Rec2 + c3)Aie−bi
ν

R2
(t+∆t) [v(t+∆t)−v(t)]

∆t

t+∆t∫

t

e
bi
ν

R2
u
du=

=(c1Re
c2+c3)Aie

−bi
ν

R2
(t+∆t) [v(t+∆t)−v(t)]

∆t

R2

biν

(
e
bi
ν

R2
(t+∆t)− ebi

ν

R2
t
)
=

=(c1Re
c2 + c3)

AiR
2

∆tbiν

(
1− e−bi

ν

R2
∆t
)
[v(t+∆t)−v(t)] (5.7)

yi(t+∆t)≈ yi(t)e−bi
ν

R2
∆t
+

(5.8)

+(c1Re
c2 + c3)

AiR
2

∆tbiν

(
1− e−bi

ν

R2
∆t
)
[v(t+∆t)−v(t)]

τwu(t+∆t)=
2µ

R

8∑

i=1

yi(t+∆t) (5.9)

τwu(t+∆t)≈
2µ

R

8∑

i=1

yi(t)e
−bi

ν

R2
∆t
+

(5.10)

+(c1Re
c2 + c3)

AiR
2

∆tbiν

(
1− e−bi

ν

R2
∆t
)
[v(t+∆t)−v(t)]

At the beginning of simulations (i.e. at the first iterative step), the value of all
terms of yi(t) equals zero. At every subsequent time step, these values change
according to Eq. (5.5).



150 Z. Zarzycki et al.

6. Results of waterhammer simulations and efficiency of the

method

6.1. Governing equations

Theunsteady flow, accompanying thewater hammer effect,may bedescri-
bed by a set of the following partial differential equations (Wylie and Streeter,
1993):

— equation of continuity
∂p

∂t
+ c20ρ0

∂v

∂x
=0 (6.1)

— equation of motion
∂v

∂t
+
1

ρ0

∂p

∂x
+
2τ

ρ0R
=0 (6.2)

where v = v(x,t) is the average velocity of liquid in the pipe cross-section,
p = p(x,t) – average pressure in the pipe cross-section, R – inner radius of
the tube, τ – wall shear sterss, ρ0 – density of fluid, c0 – velocity of pressure
wave propagation, t – time, x – axial location along the pipe.

Awaterhammer simulation is carried out by themethod of characteristics
(MOC).Thismethod (for smallMachnumber flow) transformequations (6.1),
(6.2) into the following systems of equations (Wylie and Streeter, 1993)

±dp±ρ0c0dv+
2τwa

R
dt=0 (6.3)

if

dz=±c0dt (6.4)

Details of the further procedure of calculation may be found in papers by
Zarzycki and Kudźma (2004) or Zarzycki et al. (2007).

6.2. Numerical simulation

The investigated system (tank-pipeline-valve) is presented in Fig.6.

For calculations, the following data were assumed (Wichowski, 1984):

– length of hydraulic line L=3600m

– inner diameter of tube d=0.62m

– thickness of tube wall g=0.016m

– modulus of elasticity of tube wall material (cast iron) E =1.0 ·1011Pa
– bulk modulus of liquid EL =2.07 ·109Pa



Improved method for simulating transients... 151

Fig. 6. Scheme of the pipeline system; 1 – hydraulic pipeline, 2 – cut-off valve,
3 – tank, v0 – direction of fluid flow (initial condition)

– velocity of pressure wave propagation c0 =1072m/s

– initial velocity of fluid flow v0 =1.355m/s

– density of water ρ0 =1000kg/m
3

– kinematic viscosity of fluid ν =1.31 ·10−6m2/s

The critical Reynolds number is calculated according to the following equ-
ation (Ohmi et al., 1985)

Rec =800
√
Ω (6.5)

where Ω = ωR2/ν is the dimensionless angular frequency; ω = 2π/T – an-
gular frequency; T = 4L/c0 – period of pressure pulsation, L – length of
pipe.
The value of Reynolds number at the initial condition is Re= 2Rv0/ν =

6.41·105, the critical Reynolds number according toEq. (6.5) Rec =1.48·105,
so (Re>Rec) and the flow is regarded as turbulent.
The unsteady state of the system was achieved by a complete closure of

the valve for a time longer than a half of the wave’s period in the pipe, i.e.
tz = 90s > 2L/c0 = 6.72s, consequently it can be regarded as an example
of a complex waterhammer (Wylie and Streeter, 1993). Closure of the va-
lve proceeded according to the characteristics presented in Fig.7 (Wichowski,
1984). The presented characteristics illustrate the variable liquid flow rate Qz
through the valve due to its slow shutdown.
Two simulation cases were investigated:

• in the first one, the unsteady stress τwu was calculated bymeans of Eq.
(5.10) bymaking use of new weighting function (4.6),

• in the second one, the unsteady stress τwu was calculated using the
traditional method (Zielke, 1968).

The results of numerical simulation at the valve carried out for a given
number of time steps is presented in Fig.8. As one can see, after the valve



152 Z. Zarzycki et al.

Fig. 7. Flow characteristics of closing of the water supply valve

closure which lasts 90s, there are still evident significant pressure pulsations.
This effect is caused by very large inertia of the liquid (high value of the
Reynoldsnumberandasubstantial lengthof the liquid line).Therefore,despite
the long timeof valve closure, disappearingpressurepulsations are still present
in the line.

Fig. 8. Pressure fluctuation at the valve after its closure

Thecomparison of two simulations, using the old andnewmethodof calcu-
lations of unsteady frictional losses presented inFig.8 shows a good agreement
of the results obtained, and thus confirms that the model of efficient compu-
tation satisfactorily replaces the traditional method.

However, an important difference that can be observed in the simulations
is the time of numerical calculations. Figure 9 presents the dependence of
time necessary to carry out calculations on the number of time steps used for
the simulations. As it can be clearly seen, the time of numerical calculations
carried out by means of the traditional method using Zarzycki’s weighting
function (4.1) ismuch longer than that used to calculate it using the presented



Improved method for simulating transients... 153

Fig. 9. Comparison of numerical computation times

method.Moreover, the time necessary to carry out the calculations conducted
for the traditional method increases exponentially, along with the increase of
the number of time steps.However, for the effectivemethod the time increases
slowly in a linear way.
On thebasis of the above presented investigations, it canbe stated that the

presently developed method of simulating pressure changes and flow velocity,
while taking into account the unsteady friction loss expressed by (5.10), is
an effective method which can replace the tradition method and which offers
satisfactory accuracy.

7. Conclusions

The paper presented the transfer function relating the pipe wall shear stress
and themean velocity. It is based on a four-region distribution of eddy visco-
sity at the cross-section of the pipe. This model was used to determine the
unsteady wall shear stress as integral convolution of acceleration and weigh-
ting function similarly as it was done for laminar flow (see e.g. Zielke, 1968).
However, the weighting function has a complex mathematical form and it is
not very convenient to use it for simulating transient flows. Therefore, it was
approximated to a simpler form, which was a sum of eight exponential sum-
mands. For this form of the weighting function and approximated in this way,
an effective method of simulating wall shear stress was developed.
The most important conclusions which can be drawn from this study are

as follows:

• The weighting function depends on the dimensionless time t̂ = νt/R2
(as in the case for laminar flow) and additionally on the Reynolds num-
ber Re.



154 Z. Zarzycki et al.

• Theweighting functionapproaches zeroas the time increases. For agiven
value of dimensionless time it is greater in the case of laminar flow and it
decreases along with the increase of Reynolds number Re in the case of
turbulent flow. It means that the hydraulic resistance with an increase
in the Reynolds number Re goes faster to a quasi-steady value, which is
associated with increase of damping.

• Approximatedweighting function (4.6) proves to be a precise expression
in the range 10−5 ¬ t̂ ¬ 10−1 and 2 · 103 ¬ Re ¬ 107 (relative error
smaller than 5%).

• Thepresent approximatedmethodof calculatingwall shear stress (which
incidentally is similar to Trikha’s method developed for laminar flow)
makes it possible to significantly reduce the time necessary to carry out
calculations as compared with the traditional method.

• The weighting function developed in this study shows a good fit with
Vardy’s model for short times (high frequencies). However, the discre-
pancy between the values of weighting function obtained in both me-
thods increase together with the increase of Reynolds number. It should
be pointed out that both functions were obtained on the basis of deter-
mined distribution of eddy viscosity. This is justified for high frequencies
(Brown et al., 1969). It seems that for themean band of frequencies ano-
thermodelwould bemore appropriate.Amodel inwhich eddyviscosity,
except for viscous sublayer, would be determined using turbulent kinetic
energy models (e.g. k−ε, k− l).
• The weighting function has a form which can easily undergo the first
Laplace transformation and then Fourier transform. This is useful for
determination a frequency characteristics of a long hydraulic pipeline
systems.

References

1. Abreu J.W., Almeida A.B., 2000, Pressure transients dissipative effects: a
contribution for their computational prediction, 8th International Conference
on Pressure Surges, BHRGroup, The Hague, Netherlands, 449-517

2. Bergant A., Simpson A.R., 1994, Estimating unsteady friction in transient
cavitating pipe flow,Proc., 2nd Int. Conf. On Water Pipeline Systems, BHRA
Group Conf. Series Publ., 10, UK, Edinburgh, 24-26May, 2-15



Improved method for simulating transients... 155

3. Bergant A., Simpson A.R., Vitkovsky J.P., 1999, Review of unsteady
friction models in transients pipe flow, 9th International Meeting of the Work
Group om theBehaviour ofHydraulicMachinery under SteadyOscillatoryCon-

ditions, IAHR, Brno, Czech Republic, September 1-12

4. Brown F.T., Margolis D.L., ShahR.P., 1969, Small amplitude frequency
behavior of fluid lineswith turbulent flow,Journ. of Basic Eng., Trans. ASME,
s. D, 4, 678-692

5. BrunoneB.,GoliaU.M.,GrecoM., 1994,Estimating unsteady friction in
transient cavitating pipe flow,Proceedings of the 2nd International Conference
on Water Pipeline Systems, Edinburgh 1994, BHRGroup, 1-16

6. Brunone B., Golia U.M., Greco M., 1995, Effects of two-dimensionality
onpipe transientsmodeling,Journal of Hydraulic Engineering, ASCE,121, 12,
906-912

7. Bughazem M.B., Anderson A., 2000, Investigation of an unsteady friction
model for waterhammer and column separation, 8th International Conference
on Pressure Surges, BHRGroup, The Hague, Netherlands 2000, 483-498

8. Carstens M.R., Roller J.E., 1959 Boundary-shear stress in unsteady tur-
bulent pipe flow, Journ. of Hydraulics Division, Proceedings of the American
Society of Civil Engineers, February, 67-81

9. Dailey J.W. Hankey W.L., Olive R.W., Jordan J.M., 1956, Resistan-
ce coefficient for accelerated and decelerated flows through smooth tubes and
orifices, J. Basic Eng., 78, 1071-1077

10. FerranteM., BrunoneB., 2007,Wavelets for the analysis of transient pres-
sure signals for leak detection, Journal of Hydraulic Engineering, ASCE, 133,
November, 1274-1284

11. Ghidaoui M.S., Mansour S., 2002, Efficient treatment of the Vardy-Brown
unsteady shear in pipe transients, Journal of Hydraulic Engineering, ASCE,
January, 102-112

12. Hinze J.O., 1953,Turbulence, Mc Graw-Hill, NewYork

13. Holmboe E.L., Rouleau W.T., 1967, The effect of viscous shear on tran-
sients in liquid lines, Journ. of Basic Eng., Trans. ASME, s. D,March,89, 11,
174-180

14. Hussain A.K.M.F., Reynolds W.C., 1975, Measurements in fully develo-
ped turbulent channel flow, ASME Journal of Fluids Engineering, December,
568-580

15. Kim S.H., 2008, Impulse response method for pipeline system equipped with
water hammer prediction devices, Journal of Hydraulic Engineering, ASCE,
134, July, 961-963



156 Z. Zarzycki et al.

16. Kudźma S., 2005,Modeling and simulation of dynamical runs in closed condu-
its of hydraulics systems using unsteady frictionmodel, PhDwork at Szczecin
University of Technology [in Polish]

17. Kwon H.J., Lee J.J., 2008, Computer and experimental models of transient
flow in a pipe involving backflowpreventers, Journal of Hydraulic Engineering,
ASCE, 134, April, 426-433

18. LiivU.R., 1965,OnLosses foracceleratedmotionof liquid inpressureconduits,
TrudyTallinskogo, Politkhicheskogo Instituta, Series A,223, 43-50 [inRussian]

19. Ohmi M., Iguchi M., Usui T., Minami M., 1980-1981, Flow pattern and
frictional losses in pulsating pipe flow, (Part 1-7),Bulletin of JSME, 23, 24

20. Ohmi M., Kyonen S., Usui T., 1985, Numerical analysis of transient turbu-
lent flow in a liquid line,Bulletin of JSME, 28, 239,May, 799-806

21. Pezzinga G., 1999, Quasi 2D model for unsteady flow in pipe networks, Jo-
urnal of Hydraulic Engineering, ASCE, 125, 7, 676-685

22. Rahl F.M., Barlamont J., 1996, Unsteady pipe flow simulations using
unsteady friction flow formulae, 7th International Conference on Pressure Sur-
ges and Transients in Pipeline and Open Channels, BHR Group, Harrogate,
UK, 313-322

23. SafwatH.H., Polder J., 1973,Friction-frequency dependence for oscillatory
flows in circular pipe, Journal of the Hydraulics Division, ASCE, 99, HY 11,
1933-1945

24. Satlar A.M., Chaudhry M.H., Kassen A.A., 2008, Partial blockage de-
tection in pipelines by frequency responsemethod, Journal of Hydraulic Engi-
neering, ASCE, 134, January, 76-89

25. Satlar A.M., Dicherson J.R., Chaudhry M.H., 2009, Wavelet-Galerkin
solution to the water hammer equations, Journal of Hydraulic Engineering,
ASCE, 135, April, 283-295

26. Schohl G.A., 1993, Improved approximatemethod for simulating frequency-
dependent friction in transient laminar flow, Journal of Fluids Eng., Trans.
ASME, 115, September, 420-424

27. Solodovnikov V.V., 1964, Statistical Dynamic of Linear Control Systems,
WNT,Warszawa

28. Suzuki K., Taketomi T., Sato S., 1991, Improving Zielke’s method of si-
mulating frequency-dependent friction in laminar liquid pipe flow, Journal of
Fluids Engineering, Trans. ASME, 113, December, 569-573

29. Trikha A.K., 1975, An efficient method for simulating frequency-dependent
friction in transient liquid flow, Journ. of Fluids Eng., Trans. ASME, March,
97-105



Improved method for simulating transients... 157

30. Vardy A.E., 1980, Unsteady flows: fact and friction, Proc. 3rd International
Conference on Pressure Surges, Cantenbury, UK, 1-14

31. VardyA.E,Brown J.M.B., 1995,Transient, turbulent, smoothpipe friction,
Jurnal of Hydraulic Research, 33, 4, 435-456

32. Vardy A.E., Brown J., 1996, On turbulent, unsteady, smooth-pipe friction,
Proc. of 7th International Conference on Pressure Surges, BHRA, Fluid Eng.,
289-311, Harrogate, UK, 16-18 April, 289-311

33. Vardy A.E., Brown J.M.B., 2003, Transient turbulent friction in smooth
pipe flows, Journal of Sound and Vibration, 259, 5, 1011-1036

34. VardyA.E.,BrownJ.M.B., 2004,Transient turbulent friction in fully rough
pipe flows, Journal of Sound and Vibration, 270, 233-257

35. Vardy A.E., Brown J.M.B., 2007, Approximation of turbulent wall shear
stress in highly transient pipe flows, Journal of Hydraulic Engineering, ASCE,
November, 1219-1228

36. Vardy A.E., HwangK.L., Brown J.M.B., 1993,Aweighting functionmo-
del of transient turbulent pipe friction, Jurnal of Hydraulic Research, 31, 4,
553-548

37. Wichowski R., 1984, One-dimensional theory of water-hammer phenomenon
in water-supply system conduits, Part II,Archiwum Hydrotechniki, 1/2, 39-58

38. WylieE.B., StreeterV.L., 1993,FluidTransient in Systems,PrenticeHall,
EnglewoodCliffs, NY

39. Vitkowsky J., Lambert W., Simpson A., 2000,Advances in unsteady fric-
tionmodeling in transient pipe flow, 8th International Conference on Pressure
Surges, BHRGroup, The Hague, Netherlands, 471-482

40. Wood J.D., Funk J.E., 1970, A boundary-layer theory for transient viscous
losses in turbulent flow, Journ. of Basic. Eng., Trans. ASME, December, 865-
873

41. Zarzycki Z., 1993, Friction factor in a pulsating turbulent pipe flow, Trans.
of the Institute of Fluid-Flow Machinery, 95, 179-196

42. Zarzycki Z., 1994,A hydraulic resistance’s of unsteady fluid flow in pipes, To
be published by Technical University of Szczecin, No 516, Szczecin [in Polish]

43. ZarzyckiZ., 2000,Onweighting function forwall shear stressduringunsteady
turbulent pipe flow, 8th International Conference on Pressure Surges, BHR
Group, The Hague, 529-543

44. Zarzycki Z., Kudźma S., 2004, Simulations of transient turbulent flow in
liquid lines using time-dependent frictional losses, The 9th International Con-
ference on Pressure Surges, BHRGroup, Chester, UK, 24-26March, 439-455



158 Z. Zarzycki et al.

45. Zarzycki Z., Kudzma S., Kudzma Z., Stosiak M., 2007, Simulation of
transients flow in a hydrostatic systemwith a long liquid line, Journal of The-
oretical and Applied Mechanics, 45, 4

46. ZielkeW., 1968, Frequency-dependent friction in transient pipe flow, Journal
of Basic Eng., Trans. ASME, March, 109-115

Ulepszona metoda symulacji przepływów przejściowych w przewodach

hydraulicznych

Streszczenie

Artykuł przedstawia zagadnienie modelowania i symulacji przebiegów przejścio-
wych podczas turbulentnego przepływu cieczy w przewodach ciśnieniowych. Chwilo-
we naprężenie styczne na ściance przewodu przedstawiono w postaci całki splotowej
z funkcjiwagi i przyspieszeniacieczy.Funkcjawagidlanaprężenia stycznegonaściance
przewodu zależy od czasu bezwymiarowego i liczbyReynoldsa.Ma ona zawiłą postać
matematyczną, dlatego aproksymowano ją do prostszej postaci, przydatnej do prak-
tycznych obliczeń inżynierskich. Przedstawiono efektywny sposób rozwiązania całki
splotowej, opierając się na metodzie podanej przez Trikha (1975) dla przepływu la-
minarnego. Podano zastosowanie ulepszonej metody symulacji naprężenia stycznego
dometody charakterystyk podczas uderzenia hydraulicznego.Charakteryzuje się ona
dużą efektywnością w stosunku do metody tradycyjnej.

Manuscript received May 12, 2010; accepted for print July 28, 2010