Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

53, 3, pp. 543-555, Warsaw 2015
DOI: 10.15632/jtam-pl.53.3.543

A NUMERICAL APPROACH TO THE STANDARD MODEL OF WATER

HAMMER WITH FLUID-STRUCTURE INTERACTION

Sławomir Henclik

Institute of Fluid-Flow Machinery, Polish Academy of Sciences, Gdańsk, Poland

e-mail: shen@imp.gda.pl

In the classic water hammer (WH) theory, 1D liquid flow in a quasi-rigid pipe is assumed.
When the pipe is flexible or is fixed to the foundation with elastic supports, the dynamic
fluid structure interaction (FSI) should be taken into account formore accuratemodelling of
the system behaviour. The standardmodel ofWH-FSI for a straight pipe reach is governed
by fourteen hyperbolic partial differential equations of the first order, two for 1D liquid flow
and twelve for 3Dmotion of the pipe. Thismodel is presented in the paper and an algorithm
for its numerical solution based of themethod of characteristics is proposed.Basic boundary
conditions (BC) are shortly discussed. The important condition at the junction of two sub-
pipes fixed to the foundationwith a viscoelastic support is presented in details and a general
method of its solution is proposed.

Keywords: water hammer, transient pipe flow, fluid-structure interaction, standard model,
method of characteristics

1. Introduction

Thewater hammer (WH)phenomenonhasbeenof scientific interest for over a century.Extensive
studies on this effect can be found in theworks ofWylie andStreeter (1993), Almeida andKoelle
(1992), Adamkowski (2013) orGhidaoui et al. (2005).WH is a result of a sudden change in pipe
flow conditions due to valves operation, hydraulic machinery load variation or other reasons.
During the transient the changes in liquid velocity result in pressure variations, which may
travel along the pipe producing loads of the structure (pipe, supports, hydraulic devices) and
perturbations to the system functioning or even its damage.WHcan appear in various pipe flow
systems like hydro-power installations, nuclear reactor cooling systems, engine fuel injection
pipes and others. In some cases, the pipe is stiff and does not move, so the action of the liquid
onto the pipe is considered as quasi-static.When the pipe is elastic or thewhole piping canmove
on its supports, thismotion forced by the unsteady flow, influences in reverse the flow variables,
and the dynamic fluid-structure interaction (FSI) has to be taken into account. These effects has
been examined by scientists for a few decades and the works ofWiggert et al. (1987), Tijsseling
and Lavooij (1990), Wang and Tan (1997) or Wiggert and Tijsseling (2001) can be pointed as
the reference. A specific and a complex area of potential application ofWH-FSI concepts is the
phenomenon of blood flow in human arteries, being an important area of research, discussed e.g.
by Alastruey et al. (2012).

Three main factors responsible for FSI coupling are pointed in literature. The weakest is
the friction between the pipe-wall and the liquid. Due to physics of WH-FSI, an unsteady
friction model should be used for modelling of this effect, and one of the existed models used
for classicWH, discussed e.g. byAdamkowski and Lewandowski (2006), Vitkovsky et al. (2006)
or Urbanowicz and Zarzycki (2012) could be adopted. However, the quasi-steady model is still
quite popular and often used because of its simplicity. The Poisson effect is responsible for the



544 S. Henclik

second FSI coupling factor. Pressure variations induce circumferential stresses and strains in
the pipe-wall which are accompanied by longitudinal ones. The latter propagates along the pipe
faster than the pressurewave in the liquid andproduces in the reverse action additional pressure
oscillations known as the precursor (PC) wave. The third and the strongest FSI mechanism is
the junction coupling (JC) effect which appears at pipe bends, knees, ends, valves, reductions,
etc. In fact, this phenomenon is especially important if the pipeline is able to move as a whole
structure. Due to this motion, strong interaction between the pipe and the liquid variables can
occur.Themotion of the junction canbealso responsible for coupling betweendifferentmodes of
structural vibrations of neighbouring sub-pipes. As far as the friction and Poisson couplings are
modelled by certain terms in the governing equations, JC is modelled by boundary conditions
(BC). The JCmechanism can also be responsible for producing theWHeffect by suddenmotion
of the structure.

In the classic water hammer theory a stiff and immovable pipe with one dimensional (1D)
liquid flow is assumed. FSI is considered to be a quasi-static effect and the pipe-wall elasticity
is taken into account only within the formula for the celerity of elastic waves in the liquid. Two
hyperbolic partial differential equations (PDE) of the first order are used for description of the
liquid flowwhich allows one to find the pressure p and velocity v as a function of the position x
and time t. When the dynamic FSI is taken into account, two main approaches are proposed
in the literature. The simpler consists of four PDEs, two for 1D liquid flow and the other two
for the longitudinal motion of the straight pipe reach. These four hyperbolic equations form the
four equations (4E) model ofWH-FSI. In this model, the 1D approximation neglects the radial
motion of the pipe and the liquid and is valid for low frequency assumption. When 3D motion
of the pipe is possible, additional ten equations of lateral and torsional movement produce
altogether the fourteen equations (14E) standardmodel ofWH-FSI. For numerical modelling of
the WH behaviour, the method of characteristics (MOC) is frequently used. When FSI effects
are taken into account, the finite elements method (FEM) can be alternatively exploited to the
structural equations, which produces theMOC-FEM technique. Pure FEMbased approach are
alsoproposed.Ontheotherhand, the fullyMOCbasedalgorithm isproposed in the fundamental
paper ofWiggert et al. (1987), however only general ideas of theirmethodwere presented there.
A thorough presentation of theMOCapproach is included in thework of Tijsseling (1993), who
applied his model to a 2D pipeline system. In fact, for structural analyses, the FEM technique
can be considered to bemore suitable, however using a unifiedMOC approach to the liquid as
well as to the structure is justified due to the wavy nature of the phenomenon.

In the current paper, the standardmodel ofWH-FSI is described and aMOCbasedmethod
of numerical solution for a 3D pipeline system is presented. Themain steps of the algorithm are
discussed and some details and novel approaches are especially explained. An important part of
themethod is theBCat the junction of two sub-pipes fixed to the foundationwith a viscoelastic
support, which is formulated as a differential equation of motion. The solution to this equation
is found and presented in a compact form, convenient for numerical computations. The author
does not know any other workwhere such a detailed 3D treatment of that problem is considered
and solved.

2. Assumptions and governing equations

The coordinate system and variables used in the current approach are presented in Fig. 1.

The pipe of length L, wall thickness e, diameter D is straight, prismatic, thin-walled
(e/D ≪ 1) and slender (D/L ≪ 1). The pipe material is linearly elastic and no buckling
appears. The flow velocity v is of little relativity to the elastic waves celerity c (v/c≪ 1). The
liquid is weakly compressible, linearly elastic and its density changes are small (pressures are
low relative to the bulk modulus K, p/K ≪ 1). The low frequency assumption applied means



A numerical approach to the standard model... 545

Fig. 1. Variables and coordinate system used for modelling of theWH-FSI phenomenon

neglecting of the radial inertia effects of the liquid and the pipe wall. The friction between the
pipe and the liquid is taken into account and a quasi-steady model is assumed. The structural
damping of the pipematerial and two-phase flow are assumed not to be present. Themotion of
the system for a straight pipe reach is described by four uncoupled sets of PDEs. The first one
consists of four equations for longitudinalmotion of the system.The second set of two equations
governs torsionalwaves in thepipe.Another two sets of fourPDEs each define the lateralmotion
of the pipe in two perpendicular planes – vertical (0zx) and non-vertical (0yx). It is assumed
that the liquid does not interact with the pipe in torsional vibrations. For the lateral movement,
theTimoshenko beammodel is usedwith the liquid being accounted as an addedmass for linear
motion. Liquid rotational inertia effects and the Coriolis force due to bending angular velocity
are neglected.
Four equations govern the longitudinal motion of the system. Their detailed derivation can

be found e.g. in Wang and Tan (1997). Two of the equations describe 1D liquid flow

∂v

∂t
+
1

ρ

∂p

∂x
=−g sinα−

4τs
ρD

∂v

∂x
+
1

ρc2
∂p

∂t
=2ν
∂wx
∂x

(2.1)

where ν is the Poisson coefficient and ρ is liquid density. The pressurewave celerity c is given as

c=

√
K

ρ

1
√
1+χ(1−ν2)

(2.2)

The right-side factor in the formula above expresses the influence of pipe-wall elasticity. The
parameter χ used in it is defined as

χ=
KD

Ee
(2.3)

The other two equations of the 4E model govern the longitudinal pipe motion defined with
velocity wx and stresses σx (σx =Qx/As,As is the pipe cross-section area)

∂wx
∂t
− 1
ρs

∂σx
∂x
=−gsinα+ τs

eρs

∂wx
∂x
− 1
ρsc
2
s

∂σx
∂t
=− νD
2Ee

∂p

∂t
(2.4)

The pipematerial density is ρs and the longitudinal elastic waves travel in it with the celerity cs

cs =

√
E

ρs
(2.5)



546 S. Henclik

The gravitational term (g is acceleration of gravity) depends on the angleαbetween the horizon
and the pipe. One can easily identify the components responsible for the Poisson and friction
couplings. For the current quasi-steady liquid pipe-wall friction model, the following formula is
used for determination of the shear stresses τs

τs =
λρ

8
(v−wx)|v−wx| (2.6)

In the equation above, λ is the Darcy-Weisbach friction factor and should be adequately deter-
mined. The pipe torsional vibrations are governed by the standard two equations (like for the
rod)

∂ωx
∂t
−
1

ρsI0

∂Mx
∂x
=0

∂ωx
∂x
−
1

GI0

∂Mx
∂t
=0 (2.7)

where Mx [Nm] is torsional moment of the force, ωx [1/s] – angular torsional velocity, I0 [m
4]

– pipe cross-section polar moment of inertia and G [Pa] – shear modulus of the pipe material,
G=E/[2(1+ν)]. The shape of the above equations is the result of no-interaction assumption
between the pipe and the liquid in torsional vibrations. However, if non-zero amount of accom-
panying water is going to be taken into account, then the factor ρsI0 in Eq. (2.7)1 should be
adequately increased.

In the lateral motion at the non-vertical plane the governing equations for the Timoshenko
beam model (Timoshenko and Young, 1955; Meirovitch, 1967) and the assumed coordinate
system are found to be

∂wy
∂t
+
1

m

∂Qy
∂x
=0

∂wy
∂x
+
1

T

∂Qy
∂t
=ωz

∂ωz
∂t
− 1
b

∂Mz
∂x
=−1
b
Qy

∂ωz
∂x
− 1
s

∂Mz
∂t
=0

(2.8)

The variables above are the linear velocitywy [m/s], shear forceQy [N], bending angular velocity
ωz [1/s] and the bendingmoment of force Mz [Nm]. The parameters m, b, s, T are introduced
formore clear presentation of themodel and the assumptions. They are defined in the following
way

m= ρsAs+ρAc = ρsAs(1+ ξ) T =κGAs b= ρsI s=EI (2.9)

The form of equations (2.9)1 and (2.9)3 results from the water-in-pipe accountingmethod – the
water is taken into account in the former, not in the latter. Other possibilities can be tested by
propermodelling of these parameters. InWiggert et al. (1987) the expression for b has included
the water component as well and the current assumption was proposed by Tijsseling (1993). In
fact, the influence of the weighted water component in b can be numerically tested, calibrated
and concluded. The parameter ξ in (2.9)1 expresses the ratio of the water-to-pipe mass and is
given in

ξ=
ρD

4ρse
(2.10)

In Eq. (2.9)3 I [m
4] is the axial moment of inertia of the pipe cross-section, and κ in (2.9)2 is

the shear coefficient. A simple approach gives for a circular tube the value ofκ=0.5, but due to
more detailed analyses, slightly larger values are estimated.κ=0.5+ν/(8+6ν) was proposedby
Cowper (1967) andκ=0.5+ν/(4+2ν)was estimated byHutchinson (2001). In the vertical 0xz



A numerical approach to the standard model... 547

plane, the equations are similar (the proper signs are the result of dextrose coordinate system)
and the gravity term is also included

∂wz
∂t
+
1

m

∂Qz
∂x
=−gcosα

∂wz
∂x
+
1

T

∂Qz
∂t
=−ωy

∂ωy
∂t
+
1

b

∂My
∂x
=
1

b
Qz

∂ωy
∂x
+
1

s

∂My
∂t
=0

(2.11)

3. Transformation of the equations

3.1. Longitudinal waves

The fundamental step of theMOC technique is linear transformation of the governing equ-
ations to get CE, which have the partial derivatives against x and t formed into absolute de-
rivatives with time for a certain dependence x(t). Such transformation applied to longitudinal
equations (2.1) and (2.4) results in two sets of two equations each. They govern two coupled
waves – WH and PC ones. The former propagates with the celerity c0 (equations C0) and the
equations C1 govern the PCwave which propagates with the celerity c1. Because of the Poisson
coupling, these celerities are slightly different than the original celerities c and cs given with
Eqs. (2.2) and (2.5). Moreover, these waves are not pure liquid and pure pipe waves as can be
seen below. The C0 equations are

d

dt
(v+Sw1)+ε

d

dt
(r− S̃q1)=−(1+S)g sinα− (1−R)

4τs
ρD

(3.1)

Theyarevalid for thepositiveandnegative characteristic slope (ε=±1)andthedependencex(t)
defined with dx/dt= εc0. The celerity of theWHwave is

c0 =
c√
A

(3.2)

In Eq. (3.1), unified variables are used. The unified velocity w1 is just the pipe section longitu-
dinal velocity wx. The unified pressure and longitudinal stresses are given by

r=
p

ρc0
q1 =

σx
ρsc1

(3.3)

The following pair of equations (ε=±1) referred to as C1 governs the PCwave

d

dt
(Rv−w1)+ε

d

dt
(R̃r+q1)= (1−R)gsinα− (1+S)

τs
ρse

(3.4)

The characteristic slope is dx/dt= εc1 and the PCwave celerity c1 is given as

c1 = cs
√
A (3.5)

The presented transformations are valid if the below relation holds

1−γ+χ> 0 (3.6)

where

γ =
Kρs
Eρ

(3.7)



548 S. Henclik

Relation (3.6) is valid for all practical cases, and for water in metal pipe and thin-walled as-
sumption is even stronger (> 1). It is useful to define the following parameter B

B=
1

2

[
(1−γ+χ)+

√
(1−γ+χ)2+4γχν2

]
(3.8)

The correcting parameter A, which is slightly greater than 1 (A = 1 for ν = 0) can now be
presented as

A=1+
χν2

B−χν2
(3.9)

The Poisson coupling parameters S,R are small (R=S=0 for ν =0) and defined below

S =
2γν

B
R= ξS (3.10)

The tilde S,R parameters are modified with the celerities ratio η= c0/c1 and given as

R̃=Rη S̃=
S

η
(3.11)

The developed clear form of the CE has an important advantage as it allows for preliminary
analyses of physical behaviour of the system even prior to computations. If the right sides of
them are assumed zero (no friction and a horizontal pipe) and the variables in the parenthesis
on the left sides are considered as new “wave variables” – velocities and stresses, these equations
would represent simple elastic waves. But due to the shape of the wave variables they are not
pure liquid (WH) nor pure pipe (PC) waves. Thus the real physical quantities, velocities and
stresses (pressures), are superpositions of the WH and PC waves. The dominant role is played
by the former in liquid variables and the latter in the structure ones.These effects are possible to
observe in experimental and numerical records (see Adamkowski et al., 2010 or Henclik, 2010).

3.2. Torsional vibrations

Transformation of governing equations (2.7) results in compatibility equations C4

dw4
dt
−ε
dq4
dt
=0 (3.12)

They are valid for the characteristic slope dx/dt= εc4, and the celerity c4 of torsional waves is

c4 =

√
G

ρs
(3.13)

TheCE and unified variables, velocities w and loads q, are numbered consequently with indexes
1, 2, 3 for linear and 4, 5, 6 for angular degrees of freedom, corresponding respectively to the x,
y, z axes. The variables in C4 equations are the normalized torsional angular velocity w4 and
the moment of force q4

w4 =hωx q4 =
hMx
ρsI0c4

(3.14)

The parameter h [m] introduced above to get the unified variables can be in fact arbitrarily
selected. But for the model homogeneity, its value is assumed to be the same as calculated
further for lateral motion.



A numerical approach to the standard model... 549

3.3. Lateral motion

Equations (2.8) of transversal movement in the non-vertical plane are transformed to CE
representing shearing waves C2 coupled with bending waves C6. The equations C2 are given in

dw2
dt
+ε
dq2
dt
= εΩw6 (3.15)

They are valid for the dependence x(t) having the slope equal to the wave celerity dx/dt= εc2

c2 =

√
T

Tm
(3.16)

The parameterΩ on the right side of (3.15) is a result of specific normalization of the equations
and is explained further in (3.20)2. The C6 equations have the following form

dw6
dt
−εdq6
dt
=−Ωq2 (3.17)

They are valid for dx/dt= εc6, where

c6 =

√
s

b
(3.18)

The unified variables are

w2 =wy q2 =
Qy
mc2

w6 =hωz q6 =
hMz
bc6

(3.19)

The compact form of the lateral CE has been achieved by specific selection of the constant h to
keep the same scaling coefficient Ω [1/s] on the right side of them. The result of this procedure
is

h=

√
b

m
Ω=

√
T

b
(3.20)

Estimation of expression (3.20)1 for the water in ametal pipe gives a value of aboutD/4 for h.

In the vertical plane, we get C3/C5 coupled CE. The shearing waves are governed by C3
equations

dw3
dt
+ε
dq3
dt
=−gcosα−εΩw5 (3.21)

They are valid as usual for the dependencex(t) being the path of thewave that propagates with
the same celerity as in the C2 case (c3 = c2). The bending waves C5 have the same celerity as
in the C6 case

dw5
dt
+ε
dq5
dt
=Ωq3 (3.22)

The definition of unified variables is analogous as in the C2-C6 case.



550 S. Henclik

4. Numerical method

The first step of the numerical algorithm is to integrate each CE in time within a specific time
step ∆t to get finite difference equations which can be solved for subsequent time moments
at the x-t plane. A proper construction of the numerical scheme is required. For stability and
convergence, the CFL (Courant-Fredicks-Lewy) condition (Quarteroni andValli, 1994; Almeida
and Koelle, 1992), which determines the relation between the time step ∆t and space grid
size∆x, is necessary

CN =
c∆t

∆x
¬ 1 (4.1)

It is easy to fulfill the condition for one pair of CE (one wave). KeepingCN equal to one allows
one to find the solutions without interpolation with the use of a scheme presented at the grid
in Fig. 2 for the C4 torsional wave. Integrating (3.12) within the time ∆t4, we get a set of two
linear equations

w4−q4 = b4L w4+q4 = b4R (4.2)

The “history parameters” b4 at the right side are a result of integrating the CE, and are given
in

b4L =w
(4L)
4 −q

(4L)
4 b4R =w

(4R)
4 + q

(4R)
4 (4.3)

The initial values of variables are taken at the beginning of C4 left (L) and right (R) characte-
ristics as it is presented at the right part of Fig. 2. Two kinds of average parameters are defined
and used consequently within the model, half of the sum and half of the difference

b4 =
b4R+b4L
2

◦

b4 =
b4R− b4L
2

(4.4)

Now, the solutions to (4.2) at the final point F (in Fig. 2) can be presented as

w4 = b4 q4 =
◦

b4 (4.5)

Fig. 2. A computational scheme for two coupled waves C6/C2 of the celerity ratio 3:2 (at the left) and
for a single wave C4 of the time step∆t4 = δt (at the right)

Amore complex case is for two coupledwaves that propagate with different celerities, which
is illustrated at the left part in Fig. 2 for the C2-C6 equations. If the ratio of shearing and
bendingwaves celerities is a rational number, we can integrate each pair of CE (C2/C6) within
a specific time step ∆t, (∆t2 and ∆t6 respectively), which is an integer multiplication of the
elementary step δt.



A numerical approach to the standard model... 551

In Fig. 2, the celerity ratio is 3:2 so the time step ratio is 2:3. This construction produces the
same space step ∆x for both waves and it is possible to move across the x-t plane and to solve
the equations for the unknown variables at each node without interpolation, as we can always
find the starting point in the past for any characteristics. Smallmodifications are required at the
beginningmoments, but then we can start the characteristic at any position x of the t=0 line
due to theknown initial conditions.As an alternative, a numerical schemewith interpolation can
be also considered. After integration, the following set of four linear equations for the unknown
variables is found

w2+q2−µ2w6 = b2L w2−q2−µ2w6 = b2R
w6−q6+µ6w2 = b6L w6+q6−µ6w2 = b6R

(4.6)

In the integration process, the mean values of the right-hand side variables in (3.15) and (3.17)
are assumed to be the arithmetic averages of their initial and final values, which produces the
following form of µ

µi =
Ω∆ti
2

(4.7)

The formulas for the history parameters can be easily derived and will not be presented herein.
Due to the coupling between the waves, the parameters µ are assumed to be upper bounded
for the stable numerical solutions, which produces an upper limit for the space size ∆x of the
grid mesh. An estimation made with the application of an iterative method (see Bjoerck and
Dahlquist, 1974) to partially transformed equations (4.6) has given for it a value of aboutD/2.
With slightly different analyses, Tijsseling (1993) suggested that the limit for an empty pipe
should be D/(2

√
2). In fact, the real space step will be selected with numerical tests. The

solutions to (4.6) are

w2 = b2 q2 =
−
◦

b2+µ2b6
1+µ2µ6

w6 =
b6+µ6

◦

b2
1+µ2µ6

q6 =
◦

b6 (4.8)

An analogous procedure can be applied for lateral equationsC3/C5 in the vertical plane. For
longitudinal motion, the results are calculated in a similar way, though a comment is necessary
as the friction terms on the right-hand side of Eq. (3.1) and (3.4) are non-linear. Usually, a li-
nearization scheme is proposed and a directmethod of solution is used. However, the non-linear
friction terms are small, so at the current approach they will be left within the history parame-
ters and assumed as fixed values with the possibility of iterative improvement. One additional
iteration is usually enough. Such a solution shows the advantage of using the history parameters.
It is a convenient idea as the results dodependonly on those parameters and they can bedefined
within a numerical routine in a different way for any changes in themodel. Introducingmaterial
damping, changing the pipe-wall frictionmodel or using interpolation only influences the shape
of those parameters. TheC0/C1 history parameters include also the gravitational term (C3/C5
ones as well). Finally, the solutions to C0/C1 are

v=
b0+Sb1
1+RS

w1 =
−b1+Rb0
1+RS

r=−
◦

b0+ S̃
◦

b1
1+RS

q1 =
−
◦

b1+ R̃
◦

b0
1+RS

(4.9)

Each group of equations is solved independently, but the same elementary time step δt is
necessary for all of them to keep the effectiveness of the method. This is because the coupling
between variables of different groups may occur at boundaries, thus common time nodes are
required there. The common space nodes at the boundaries can be also useful. In this case,
the ratios of all the wave celerities should be rational numbers. Such a condition is always true



552 S. Henclik

within a certain precision, but the algorithm is effective if the integers for these ratios are not
too large, so a wave-speed adjustment can be used by slightly changing the selected material
or geometrical parameters of the system. This procedure is also used in other approaches, and
the current author has developed and tested an algorithm allowing one to adjust all the five
celerities of the elastic waves.

5. Boundary conditions

The piping is modelled as a set of straight pipes connected at junctions where certain relations
are valid as the boundary conditions. In a more general case, some working elements may be
also taken into account. The termnon-pipe element (NPE) is used inAlmeida andKoelle (1992)
for them and theymay be a valve, pump, turbine or other hydraulic device. NPE is assumed to
be a rigid andmassive element and can be fixed to the foundationwith rigid or elastic supports
(or staying unfixed). It is connected to individual pipes by a number of inlets-outlets for which
hydraulic characteristics are determined.Asimple example ofNPE is avalvewith its dependence
between the pressure at the inlet at the left (pL) and the outlet at the right (pR) given as

pL = pR+
1

2
ρ(v−wx)2ζ(t) (5.1)

In the equation above, the loss factor ζ(t) is dependent on the valve opening degree, so it may
change in time. For a complete valve closure, it becomes infinity what in fact changes the BC
to have the form v = wx. NPE will not be considered within the current model, however the
junction is in fact a simple NPE of a constant and small loss factor. The mass and size of the
junction canusually beneglected, but the general case canbe analyzed aswell.Using the current
assumptions, the main dynamic liquid BC at the junction of two sub-pipes defines the pressure
balancewith a relation similar to (5.1). The condition is even simpler as the junction losses ζ are
constant and usually small, which produces a possibility for iterative solution. A fundamental
BC for the liquid is the continuity equation. For the FSI case and two sub-pipes at the junction,
it has the following form

v(L)−w(L)x = v(R)−w(R)x (5.2)

In spite of the fact that the junction is rigid, the velocities of the left and right edge of it used
above are not, in general, the same as 3D motion is considered and the left and right pipe
coordinate systems can be rotated (when the junction is a bend). The rigidity of the junction
means it is rigid in itself and is rigidly fixed to the pipes. That is why the junction velocity w
and pipe ends velocities w(L) and w(R) are uniquely related. In fact, they are identical for the
dimensionless junction

w
(L) =w(R) =w (5.3)

In this case, w is 6D velocity with 4,5,6 components being angular velocities multiplied by h
defined in (3.20)1.Theabove isavector equation, sowhenusing its representation inacoordinate
system, proper transform matrices of directional cosines U(L), U(R) have to be applied. If the
junction size is not negligible, the above condition has to be modified because linear motion of
the junction edge (the pipe end) has a component being the result of the junction rotation as a
rigid body.

Themotion of the junction is dependent on the way it is fixed to the foundation. For a rigid
fixing, there is no motion and the BC is w= 0. When the junction is unfixed or fixed with an



A numerical approach to the standard model... 553

Fig. 3. Junction of two pipes fixed with a viscoelastic support

elastic support as it is presented in Fig. 3, the following junction equation of motion (EOM)
should be formulated

Mü+Cu̇+Ku=FL+FP (5.4)

This is 6D equation for the junction displacement u having three linear (1,2,3) and three
angular (4,5,6) coordinates. At the right-hand side of the EOM, there are pipe and liquid forces
(andmoments) acting onto the junction. On the left-hand side, there is an inertia termwith the
matrixM and the forces from the elastic supportwith stiffness and dampingmatricesK andC.
The Kelvin-Voigt model for the support reaction is assumed. In the EOM, the convention of
multiplying angular (4,5,6) coordinates in the vector u by the parameter h is still valid and
another one for normalization of the angular coordinates of the forces (themoments) bydividing
them by h is used as well. An adequate normalization is also done to thematrices M,K,C on
the left-hand side. In general, determining these matrices is a separate task, but this problem
will not be discussed herein. The EOM is simplified in special cases, e.g. for massless junction
M=0, for zero dampingC=0 and for not supported junctionK=0,C=0.
The forces on the right side of the EOMdepend on the flow and pipemotion variables. The

liquid forces FL are mainly the result of pipe static pressures and the pipe forces FP appear
explicitly in theCE, though proper transformations are required. Form (5.4) of the EOMallows
one to understand the junction couplingmechanism, however it is not convenient for numerical
solution. The right-hand side of the EOM can be transformed using the difference form of CE
discussed in Section 4 to express the liquid and pipe forces as a linear function of the junction
velocity w= du/dt in a new time instant. Finally, the EOM can be presented in the following
form

M̃ü+ C̃u̇+K̃u=−Pw+a (5.5)

In the above equation, the matrices on the left-hand side (with tilde) are the original ones
dividedby (ρsAscs).The essence of themethod is the couplingmatrixPwhich is symmetric and
positive definite. It depends on the pipe-junction geometry and liquid-pipematerial parameters.
The junction history vector a depends on the previous values of the system variables at the
junction and its neighborhood, and changes in each time step. Both formulas for the matrix P
and vector a have been determined by the author to be applied in the algorithm. To find the
solution toEq. (5.5), its left-hand sidehas to be transformed to afinite difference form.Applying
theNewmarkmethod to it allows one to find the junction velocityw in a new time instant (the
variables with “0” indexes are taken at the previous time instant)

w=
( 2
∆t
M̃+P+ C̃+

∆t

2
K̃
)
−1[
a+M̃

( 2
∆t
u̇0+ ü0

)
−K̃
(
u0+

∆t

2
u̇0

)]
(5.6)

Thematrix in the brackets is constant for the junction and positive definite, thus the equation
can be easily solved in each time step. Knowing the velocityw, all other variables can be found.



554 S. Henclik

6. Summary

The water hammer phenomenon with dynamic fluid-structure interaction is discussed in the
paper. The standard mathematical model of WH-FSI and an algorithm for its numerical so-
lution on the basis of the MOC technique is presented. The fourteen governing equations are
transformed into compatibility equations developed in a convenient, compact formwith the use
of “unified variables” – velocities and loads. Themain steps and ideas of the numerical method
are discussed. The concept of using “history parameters” allows one to effectively design the
algorithm. The numerical scheme uses a properly designed computational grid exploiting the
technique of wave-speed adjustment but the scheme with interpolation is also possible. An im-
portant part of the algorithm is the boundary condition at the junction fixed to the foundation
with a viscoelastic support. It is formulated as a differential equation of motion, transformed
and effectively solved with the use of a special coupling matrix and the Newmark method. A
similar approach has been applied to the 1D 4E model and it works correctly. The algorithms
are now being implemented in a computer code, and the numerical results will be compared
with the records from experiments which are planned for the real 3D pipeline system built in
the laboratory.

Acknowledgements

The results presented in the paper have been partially developedwithin research project No. NN504

478839 funded by theMinistry of Science and Higher Education of Poland.

References

1. Adamkowski A., 2013,Liquid Unsteady Flows in Closed Conduits (in Polish), Publishing House
of the Institute of Fluid-FlowMachinery, Gdansk

2. Adamkowski A., Henclik S., Lewandowski M., 2010, Experimental and numerical results of
the influence of dynamic Poisson effect on transient pipe flow parameters, 25th IAHR Symposium
on Hydraulic Machinery and Systems, Timisoara, Institute of Physics Conference Series: Earth

and Environmental Sciencees, DOI: 10.1088/1755-1315/12/1/012041

3. AdamkowskiA.,LewandowskiM., 2006,Experimental examinationofunsteady frictionmodels
for transient pipe flow simulation,ASME Journal of Fluids Engineering, 128, 11

4. Alastruey J., Parker K.H., Sherwin S.J., 2012, Arterial pulse wave haemodynamics,Proce-
edings of the 11th International Conference on Pressure Surges, Lisbon, 401-442

5. Almeida A.B., Koelle E., 1992,Fluid Transients in Pipe Networks, ComputationalMechanics
Publications, Boston, London

6. Bjoerck A., Dahlquist G., 1974,Numerical Methods, Prentice Hall, chap. 6.9.1 (Polish trans-
lation, 1987, PWN,Warsaw)

7. Cowper G.R., 1966, The shear coefficient in Timoshenko’s beam theory, Journal of Applied Me-
chanics, 33, 6, 335-340

8. GhidaouiM.,MingZao,McilnnisD.,AxworthyD., 2005,A reviewofwater hammer theory
and practise,Applied Mechanics Reviews, 58, 1, 49-76

9. Henclik S., 2010,Mathematical model and numerical computations of transient pipe flows with
fluid-structure interaction,Transactions of the Institute of Fluid-Flow Machinery, 122, 77-94

10. Hutchinson J.R., 2001, Shear coefficients for Timoshenko beam theory, Journal of Applied Me-
chanics, 68, 1, 87-92

11. Meirovitch L., 1967,Analytical Methods in Vibrations, TheMacmillan Company, NewYork



A numerical approach to the standard model... 555

12. Quarteroni A., Valli A., 1994, Numerical Approximation of Partial Differential Equations,
Springer-Verlag, Berlin

13. Tijsseling A.S., Lavooij C.S.W., 1990,Waterhammer with fluid-structure interaction,Applied
Scientific Research, 47, 273-285

14. Tijsseling A.S., 1993, Fluid – structure interaction in case of waterhammer with cavita-
tion, PhD thesis, Delft University of Technology, Report 93-6, Faculty of Civil Engineering,
http://www.win.tue.nl/∼atijssel/pdf files/Tijsseling 1993.pdf

15. Timoshenko S., Young D.H., 1955,Vibration Problems in Engineering, Nostrad Co. Inc., NY

16. UrbanowiczK., Zarzycki Z., 2012,New efficient approximation ofweighting function for simu-
lations of unsteady friction losses in liquid pipe flow,Journal of Theoretical andAppliedMechanics,
50, 2, 487-508

17. Vitkovsky J.P., Bergant A., Simpson A.R., Lambert M.F., 2006, Systematic evaluation of
one-dimensional unsteady friction models in simple pipelines, Journal of Hydraulic Engineering,
132, 7, 696-708

18. WangZ.M.,TanS.K., 1997,Coupledanalysis of fluid transient and structural dynamic responses
of a pipeline system, Journal of Hydraulic Research, 35, 1, 119-131

19. Wiggert D.C., Hatfield F.J., Stuckenbruck S., 1987, Analysis of liquid and structural
transients in piping by the method of characteristics,ASME Journal of Fluids Engineering, 109,
161-165

20. Wiggert D., Tijsseling A., 2001, Fluid transients and fluid-structure interaction in flexible
liquid-filled piping,Applied Mechanics Reviews, 54, 455-481

21. Wylie E.B., Streeter V.L., 1993,Fluid Transients in Systems, Prentice-Hall, NJ

Manuscript received August 8, 2014; accepted for print November 4, 2014