Energy dissipation in unsteady turbulent pipe flows caused by water hammer

Computers & Fluids 73 (2013) 124–133

Contents lists available at SciVerse ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Energy dissipation in unsteady turbulent pipe flows caused by water hammer A. Riasi a,⇑, A. Nourbakhsh b, M. Raisee b a b

Department of Mechanical Engineering, University of Tehran, P.O. Box 11365-4563, Tehran, Iran Department of Mechanical Engineering, University of Tehran, Tehran, Iran

a r t i c l e

i n f o

Article history: Received 3 August 2009 Received in revised form 24 October 2012 Accepted 14 December 2012 Available online 8 January 2013 Keywords: k–x Turbulence model Water hammer Energy dissipation Dimensional analysis Turbulent production

a b s t r a c t Energy dissipation and turbulent kinetic energy production and its dissipation in unsteady turbulent pipe flows due to water hammer phenomena are numerically studied. For this purpose, the two-dimensional governing equations of water hammer are solved using the method of characteristics. A k–x turbulence model which is accurate for two-dimensional boundary layers under adverse and favorable pressure gradients is applied. The numerical results are in good agreement with the experimental data. Through an order of magnitude analysis, two dimensionless parameters have been identified which can be used for the evaluation of viscous and turbulent shear stress terms. The influence of these non-dimensional parameters on pressure oscillations, wall-shear-stress, dissipation rate as well as profiles of velocity, turbulent production and dissipation are investigated. The non-dimensional parameter P, which represents time scale ratio of turbulence diffusion in the radial direction to the pressure wave speed, is used to study the structure and strength of turbulence. It is found that for the case of P  1, for which the values of the non-dimensional groups are larger, the peaks of turbulence energy production and dissipation move rapidly away from the wall and turbulence structure is significantly changed. For the case of P  1, for which the values of non-dimensional parameters are smaller, these variations are found to be small. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction In pressurized pipes and channels, perturbation in the fluid flow is produced due to several operations such as starting or failure in pump and turbine and also fast opening or closing of the valve in order to regulate the flow rate. These actions result in a series of positive and negative pressure waves which move along the passage with wave speed. Regarding these positive and negative pressure waves, several problems such as pipe rupture, hydraulic equipment damage and cavitations might occur. Water hammer problems are usually simulated using one-dimensional water hammer equations based on the quasi-steady friction model. The main assumption of this method is that the head loss during water hammer phenomenon is equal to the head loss obtained from the steady flow. However, this assumption is not valid for most of the water hammer problems due to existence of strong gradients and reverse flows near the pipe wall. The main focus of the present paper is energy dissipation in water hammer flows. Energy dissipation in transient pipe flow was studied by Silva-Araya and Choudhry [1]. They defined the energy dissipation factor as a ratio between the energy dissipation at any instant and the energy dissipation obtained by assuming ⇑ Corresponding author. E-mail address: [email protected] (A. Riasi). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2012.12.015

quasi-steady conditions. They used this factor to improve the accuracy of one-dimensional equations of water hammer. The attenuation in laminar water hammer problems was numerically studied by Wahba [2]. He introduced a non-dimensional parameter which controls the viscous shear stress term. Moreover, he investigated the influence of this parameter on pressure damping, instantaneous wall-shear-stress and Richardson annular effect. To gain deeper insight of energy dissipation and structural variation of turbulence in turbulent water hammer problems, it is necessary to employ two-dimensional equations for water hammer simulation using accurate turbulence modeling. For this purpose, Vardy and Hwang [3], Silva-Araya and Choudhry [1,4], Pezzinga [5], Zhao and Ghidaoui [6] and Wahba [7] used several algebraic turbulence models. Recently, Zhao and Ghidaoui [8] used lowReynolds number k-e model of Fan et al. [9] for unsteady turbulent pipe flows. More recently, Riasi et al. [10] used low-Reynolds k–x turbulence model to simulate unsteady turbulent pipe flows emerging from water hammer problems. The main objectives of the present work are: (i) to introduce the non-dimensional groups for controlling viscous and turbulent shear stress terms, (ii) to investigate the influence of these nondimensional groups on pressure oscillations, wall-shear-stress, dissipation rate and velocity profiles and (iii) to study the radial distribution of turbulent kinetic energy production and its dissipation in one wave cycle of water hammer.

A. Riasi et al. / Computers & Fluids 73 (2013) 124–133

125

Notations a Cl D e E f Ftur Fvis Fw g gz GCI h H H0 k L M MV Nr Nx p P Q r Re t Td TW

acoustic wave velocity constant in eddy-viscosity pipe diameter error in non-dimension form validation comparison error Darcy–Weisbach friction factor turbulent force viscous force water hammer force gravity acceleration geopotential energy per unit mass grid convergence index specific enthalpy piezometric head initial piezometric head at valve kinetic energy of turbulence fluctuation per unit mass pipe length Mach number Mach number based on the radial velocity number of computational grids in radial direction number of computational grids in longitudinal direction observed order of method turbulence production in Eq. (27), non-dimensional parameter in Eq. (28) energy transfer by heat radial coordinate measured from the pipe center line Reynolds number Time diffusion time scale water hammer wave time scale

2. Mathematical modeling and numerical procedure As shown by Vardy and Hwang [3], the continuity and momentum equations for transient flow in pipes are written as

@H @H a2 @u a2 1 @r t þu þ þ ¼0 @t @x g @x g r @r

ð1Þ

   @u @u @u @H 1 @ @u þu þt ¼ g þ r qm  qu0 t0 @t @x @r @x rq @r @r

ð2Þ

where t = time; x = distance along the pipe; H(x, t) = piezometric head; u(x, r, t) = axial velocity; t(x, r, t) = radial velocity; u0 = stream wise velocity fluctuation; t0 = cross stream velocity fluctuation; a = acoustic wave velocity; g = gravity acceleration; q = density and m = kinematic viscosity. Here the piezometric head is defined as H(x, t) = p(x, t)/c + z, in which the p(x, t) = the average of pressure in pipe, c = qg = specific weight of fluid and z = pipe level. The main assumptions of Vardy and Hwang [3] in the derivation of these equations are: (a) the pipe wall is approximately inextensible, and thus there is no radial velocity at the pipe wall and (b) the pressure at any cross-section is uniform. Moreover, the nonlinear convective term in Eq. (2) is neglected due to small Mach number. This point will be clearly explained in Section 3. According to the Boussinesq approximation, the turbulent shear stress is approximated as

qu0 t0 ffi qmt

@u @r

where mt = eddy viscosity.

u u0 us unum U

t t0 V V0 Vol W x y z

longitudinal velocity in Eq. (1), internal energy in Eq. (19) fluctuation velocity in longitudinal direction friction velocity standard uncertainty in the simulation solution longitudinal velocity scale radial velocity fluctuation velocity in radial direction velocity, radial velocity scale initial bulk velocity volume energy transfer by work longitudinal coordinate distance measured from the pipe wall (R–r) pipe level

Greek symbols ax, bx, bk, rx, rk closure coefficients for k-x model a, b non-dimensional groups e dissipation per unit mass U dissipation function c specific weight of the fluid m kinematics viscosity mt eddy viscosity q density s shear stress sW wall shear stress sWss steady state wall shear stress X specific dissipation rate x1, x2 weighting coefficients n real positive parameter

Eqs. (1) and (2) form a system of hyperbolic–parabolic partial differential equations. Vardy and Hwang [3] employed a hybrid solution scheme in which the hyperbolic part and the parabolic part of the governing equations are respectively solved using the method of characteristics and finite difference. With linear combination of Eqs. (1) and (2) the following characteristic equations, describing perturbations movement at speeds +a and a relative to the fluid and along the positive and negative characteristics (see Fig. 1), are formed [3]:

Cþ :

dH a du a2 1 @ðrtÞ a 1 @ðr sÞ dx þ ¼ þ along ¼ þa dt g dt g r q @r dt g r @r

ð3Þ Fig. 1. Positive and negative characteristics for water hammer problem.

ð4aÞ

126

A. Riasi et al. / Computers & Fluids 73 (2013) 124–133

C :

dH a du a2 1 @ðr tÞ a 1 @ðr sÞ dx  ¼  along ¼ a dt g dt g r q @r dt g r @r

where total shear stress express as

ð4bÞ

s ¼ qm @u  qu0 t0 . @r

In unsteady turbulent pipe flows, the use of a high resolution grid is essential for accurate calculation of important variables such as velocity, turbulent kinetic energy and specific dissipation rate. Therefore, as shown in Fig. 2, a non-uniform mesh consisting fine enough nodes in the near-wall region is used. To this end the clustering factor is considered to be 1.05. The number of grid nodes in the radial and longitudinal directions is Nr and Nx, respectively (see Fig. 2). Regarding the characteristics method (Fig. 1), the Courant number, defines as C = aDt/Dx, is 1.0 and thus the time step is obtained from Dt = Dx/a. According to Fig. 1, the characteristic equations are integrated along the positive and the negative characteristics. The terms in the right hand side of Eq. (4) are approximated as weighted average of the relevant values at successive times n and n + 1 along the characteristic lines and thus, the following equations are obtained [6] nþ1 C þ : Hnþ1  x1 C t1 ðjÞðr tÞnþ1  x2 C u1 ðjÞunþ1 i;j1 þ x1 C t2 ðjÞðr tÞi;j i i;j1   a n nþ1 nþ1 þ þ x2 C u2 ðjÞ ui;j  x2 C u3 ðjÞui;jþ1 ¼ Hi1 g

Þni1;j1

þ ð1  x1 ÞC t1 ðjÞðr t

 ð1  x1 ÞC t2 ðjÞðrt   a þ ð1  x2 ÞC u1 ðjÞuni1;j1 þ  ð1  x2 ÞC u2 ðjÞ uni1;j g ð5aÞ

nþ1  x1 C t1 ðjÞðr tÞnþ1 C  : Hnþ1 i i;j1 þ x1 C t2 ðjÞðr tÞi;j   a þ x2 C u2 ðjÞ unþ1 þ x2 C u1 ðjÞunþ1 þ x2 C u3 ðjÞunþ1 i;j1  i;j i;jþ1 g

¼ Hniþ1 þ ð1  x1 ÞC t1 ðjÞðrtÞniþ1;j1  ð1  x1 ÞC t2 ðjÞðr tÞniþ1;j   a  ð1  x2 ÞC u2 ðjÞ uniþ1;j  ð1  x2 ÞC u1 ðjÞuniþ1;j1  g  ð1  x2 ÞC u3 ðjÞuniþ1;jþ1

ð6Þ

C u2 ðjÞ ¼ C u1 ðjÞ þ C u3 ðjÞ rj ¼

j X

Dr m

m¼1

r j ¼ ðr j1 þ rj Þ=2 where Drm = rm–rm1. At a given time t and a location x along the pipe, Eqs. (5a) and (5b) are valid for each rj and consequently, 2Nr linear simultaneous nþ1 characteristics provided for 2Nr unknowns Hnþ1 ; unþ1 i i;1 ; ui;2 ; . . . ; nþ1 nþ1 nþ1 nþ1 ui;Nr ; ti;1 ; ti;2 ; . . . ; ti;Nr 1 . A linear system of equations AX = B is formed, where A2Nr2Nr = coefficient matrix, X2Nr1 = unknowns vector and B2Nr1 = right hand side vector. Since the method of Vardy and Hwang [3] is very time consuming, Zhao and Ghidaoui [6] converted the linear system of equations AX = B into the two sub linear systems of equations with tridiagonal coefficient matrix. 2.1. Boundary conditions

Þni1;j

þ ð1  x2 ÞC u3 ðjÞuni1;jþ1

a2 Dt 1 g r j Drj aDtðm þ mtj1 Þ 1 rj1 C u1 ðjÞ ¼ r j Drj ðrj  rj1 Þ g aDtðm þ mtj Þ 1 rj C u3 ðjÞ ¼ rj Dr j ðr jþ1  rj Þ g C t1 ðjÞ ¼ C t2 ðjÞ ¼

ð5bÞ

where x1 and x2 are the weighting coefficients, superscript n is the temporal location and subscript i and j are indexes for spatial location. Coefficients Ct1, Ct2, Cu1, Cu2 and Cu3 are defined as follows

At the inlet, the head is predetermined and the radial velocity component is set to zero. For axial velocity components, Nr negative characteristics (Eq. (5b)) are valid and must be simultaneously solved. At the closed valve, all axial velocity components are zero and Nr positive characteristics (Eq. (5a)) are solved simultaneously to nþ1 nþ1 nþ1 yield Hnþ1 N ; tN;1 ; tN;2 ; . . . ; tN;Nr 1 . In the case of the valve is not suddenly closed, the linear reduction in time, has been considered for the axial velocity components. This reduction starts from the steady state velocity to the zero velocity. It is worth mentioning that in this case, the valve is closed very fast but not suddenly. To illustrate the implementation of the wall boundary condition, the no-slip boundary condition is imposed at the solid surface aligned along rNr = R. Thus,

ti;Nr ¼ 0:0 ui;Nr þ1 ¼ ui;Nr

ð7Þ

At the pipe axis, zero gradient boundary condition is employed for the axial velocity component and the radial velocity component is set to be zero.

Fig. 2. Introduction of variables in a part of pipe for implicit characteristics method.

127

A. Riasi et al. / Computers & Fluids 73 (2013) 124–133

Turbulent dissipation is the rate at which turbulence kinetic energy is converted into thermal internal energy.

2.2. Stability condition Using Von Neumann analysis [11], the following condition should be satisfied for the stability of numerical scheme [6]

0:5 6 x2 6 1; 0:5 6 x1 6 1

ð8Þ

or

0:0 6 x2 6 0:5; Dt 6

1 Dr 2 ð2  4x2 Þ ðm þ mt Þ

ð9Þ

2.3. Turbulence modeling The Wilcox’s k–x model [12] is employed for modeling of turbulence. It is well known that this turbulence model is fairly accurate for 2D boundary layer under pressure-gradient [12]. Furthermore, the k–x model can be integrated through the viscous sublayer without using any viscous damping. The k and x equations and eddy-viscosity formulation in Wilcox’s k–x model are written as

   2 @k 1 @ @k @u þ mt ¼ rðm þ rk mt Þ  bk kx @t r @r @r @r     @x 1 @ @x x @u 2  bx x2 rðm þ rx mt Þ ¼ þ ax mt r @r @t @r k @r k mt ¼

ð10Þ

x

where k = kinetic energy of turbulence fluctuation per unit mass and x = specific dissipation rate. The following closure coefficients are provided by Wilcox [12]

5 9

ax ¼ ; bx ¼

3 9 1 1 ;b ¼ ; rx ¼ ; rk ¼ 40 k 100 2 2

ð11Þ

3. Non-dimensional groups To obtain the non-dimensional form of turbulent water hammer equations, the following normalized variables are considered

u  t  x  r t q ; t ¼ ; x ¼ ; r ¼ ; t ¼ ; q ¼ ; q0 U L D nL=a V H u0 t0 ; u0 ¼ ; t0 ¼ H ¼ aU=ng us us

u ¼

ð13Þ

where the superscript  is used to denote dimensionless quantities, n = a real positive parameter equal to 4 for a full wave cycle (i.e., T = 4L/a) and for rapid transient flow set to one or lower value, q0 = density of the fluid in undisturbed state, U = axial velocity scale and V = radial velocity scale. Since us = (fU2/8)1/2, the order magnitude of friction velocity is (0.1U). Note that here the Joukowsky head (aU/ng) is used for the normalization of pressure head, whereas the pressure in high Reynold number flows are normalized by qU2/2 and in low Reynold number flows (creeping flows) by qmU/D. Different time scales can be identified for water hammer phenomenon as follows: TW = L/a = water hammer wave time scale and Td = D2/m = diffusion time scale. Here, the water hammer time scale is used for time normalization. Applying the above scales to the momentum equation (Eq. (2)), yields         @u @H L M 1 @  @u  @u   @u þ ðnM Þ t ¼  þ n r q V  þ ðnMÞu D Re r  q @r  @t @x @r  @x @r |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Viscousshearstressterm    L 1 @  þ 0:01n M    q u0 t0 Þ ð14Þ D r q @r |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Turbulentshearstressterm

In typical water hammer problems, as will be shown in the subsequence section, the Mach number is very small and thus convective terms can be neglected in the k and x equations. To solve k and x equations, the diffusion and the source terms are discretized using the implicit central difference expression and the time derivative is discretized using the first order forward differencing scheme. Finally, two tridiagonal matrix systems are solved for the k and x in each time step. At the pipe axis, symmetric boundary conditions for k and x are assumed: ok/or = 0, ox/or = 0, and for wall boundary condition, the kinetic energy of turbulence, k, is set to zero. In order to facilitate integration of the k–x turbulence model through the viscous sublayer, one can show that the specific dissipation rate, x, must satisfy the following asymptotic solution close to the wall: x ? 6m/ by2. In all computation, in order to resolve the variation of important variables within the viscous sublayer; the y+ = usy/m value of the first grid node adjacent to the wall is less than 1.0.

where M = U/a = Mach number, MV = V/a = modified Mach number based on the radial velocity, Re = UD/m = Reynolds number. Vardy and Hwang [3] showed that the radial velocity is of order [104  103] time the axial velocity. For most industrial applications, the orders of magnitude of the above non-dimensional groups are

L ¼ Oð102  104 Þ; M ¼ Oð104  102 Þ; D MV ¼ Oð108  105 Þ; Re ¼ Oð102  106 Þ

The change in density in unsteady compressible flows is of the order of the Mach number [13]. Since the Mach number is much smaller than unity for water hammer problems, it can be concluded that q  1. According to Eq. (15), the order of magnitude for the terms of Eq. (14) can be written as @u þ @t

ðnMÞ |ffl{zffl}

Oð104 102 Þ

2.4. Turbulent production and dissipation To study the structure and behavior of turbulence, the production and dissipation of the turbulent kinetic energy are defined as [12]

 2 @u P ¼ mt @r

ð12Þ

e ¼ bk kx The turbulent production represents the rate at which kinetic energy is transferred from the mean field to the fluctuation field.

ð15Þ

u

       @u L @u @H L M 1 @ @u þ n M V t  ¼   þ n r  q  @x @r @x @r D D Re r  q @r  |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} Oð106 101 Þ

Oð108 1Þ

  L 1 @ þ 0:01n M    ðq u0 t0 ÞÞ D r q @r |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl}

ð16Þ

Oð104 1Þ

From Eq. (16), it can be concluded that: (1) the convective terms can be neglected because of small Mach number, (2) the viscous and the turbulent shear stress terms are controlled by the non-dimensional groups nLM/DRe and 0.01nLM/D, respectively, (3) these terms are much stronger in longer pipe of smaller diameter and also in unsteady pipe flows with higher Mach number and (4) the viscous and turbulent shear stresses are more significant for larger values of n (i.e., n > 4).

128

A. Riasi et al. / Computers & Fluids 73 (2013) 124–133

Note that higher Mach numbers occur in pipes with low elasticity modules like polyethylene pipes where the pressure wave speed is less than the metallic pipes. Hereafter, these non-dimensional groups are defined as

c:v :

ð17Þ

Different interpretations can be presented for the non-dimensional groups a. For example, it can be interpreted as ratio of the different time scales:

a¼n

Lm aD2

¼

TW Td

ð18Þ

Eq. (18) shows that the viscous non-dimensional group, a, is defined as ratio of the water hammer wave time scale to the diffusion time scale. Accordingly, a and b can be represented as ratio of different forces as:

qðmU=DÞðpLDÞ F v is ¼ qðaU=nÞðpD2 Þ F W LU qu0 v 0 ðpLDÞ F tur ¼ b ¼ 0:01n ¼ Da qðaU=nÞðpD2 Þ F W a¼n

Lm

aD2

¼

ð19Þ

where Fvis = viscous force, Fw = water hammer force and Ftur = turbulent force. The reference surface for the viscous and the turbulent forces is the internal surface of the pipe while the reference surface for the water hammer force is the cross section area of the pipe. Note that in the above equation, it is assumed that u0  0.1U and t0  0.1U. These non-dimensional groups are consistent with those already reported by Wahba [2] for the non-dimensional parameter of laminar water hammer problems.

4. Energy dissipation To study the energy dissipation in water hammer problems, the energy balance should be considered for the standard test rig of water hammer including reservoir, pipe and valve. For this purpose, a control volume is applied to the pipe (Fig. 3). Accordingly, the balance of energy for this control volume is d dt

Z c:v :

! Z V2 _  uþ þ gz qdVol ¼ Q_  W 2 c:s:

V2 hþ þ gz 2

!

q~ V  d~ A

ð20Þ

where u = internal energy, gz = geopotential energy per unit mass, _ account for the net rate of energy h = the specific enthalpy and Q_ ; W transfer by heat and work across the boundary of control volume. Since the pipe wall is almost adiabatic, the heat loss is zero and this _ represents the term can be neglected from the energy equation. W work done on or by the control volume due to boundary displacement. The displacement happens because of pipe expansion and contraction during traveling of high and low pressure waves along the pipe. By integrating the energy equation during the simulation time and neglecting the gravity term, one obtains

 Z

Z c:v :

t¼T

¼ Wjt¼0  Wjt¼T 

Lm

L M a¼n 2¼n D Re aD LU L b ¼ 0:01n ¼ 0:01n M Da D

! V2 uþ qdVol 2

Z

T

"Z

0

c:s:

! V2 uþ qdVol 2 t¼0 # ! 2 V V  d~ A dt hþ q~ 2

ð21Þ

where T = computational time. In one cycle of water hammer and during time interval [0, L/a], the high pressure wave moves upstream, bringing the fluid to rest, compressing it and expanding the pipe. When the wave reaches the upstream end of the pipe (t = L/a), all the fluid is under additional head (DH), all the momentum has been lost and all kinetic energy has been converted into elastic energy. In this period, the flow direction is from reservoir to the valve and some work is done by the fluid on the boundary of the control. During time interval [L/ a, 2L/a], the flow direction is changed and directed from the valve to the reservoir. The pressure returns to the value which was normal before closure and the pipe wall returns to its normal form. In this period some work is done on the control volume with boundary displacement. For the perfect elasticity, the value of this work is approximately equal to the work transferred in time interval [0, L/a]. A similar scenario is repeated for the traveling of low pressure wave during time interval [2L/a, 3L/a] and [3L/a, 4L/a]. Therefore, the net work transfer during a cycle of water hammer is almost zero and also the net discharge flow from reservoir to the pipe is approximately zero. As a result, the right hand side of Eq. (21) is almost zero when the simulation time, T, is multiple of the wave period. Rearranging Eq. (21), the energy equation can be expressed as

Z c:v :

¼

uqdVol Z c:v :

 t¼T

Z c:v :

! V2 qdVol 2

uqdVol

t¼0



Z

t¼0

c:v :

! V2 qdVol 2

ð22Þ

t¼T

Note that the simulation time, T, is multiple of wave period. From Eq. (22), it is clearly seen that the change in the internal energy of the pipe flow is balanced with the variation of the kinetic energy of the pipe flow. The internal energy of the pipe flow is changed because of the viscous and turbulent dissipation. The dissipation function for two-dimensional turbulent-boundary-layer flows is approximated as [14]

/¼

  @u @u qm  qu0 t0 @y @y

ð23Þ

where, y = R–r, Using Boussinesq approximation, Eq. (23) can be expressed as

 2  2 @u @u / ¼ qm þ qmt @y @y |fflfflfflfflfflffl{zfflfflfflfflfflffl} |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} Viscousdissipation

ð24Þ

Turbulentdissipation

Eq. (24) represents the mean energy loss per unit volume and per unit time. Also, in Eq. (24), the turbulence dissipation has been approximated by the turbulence production. Integrating Eq. (24) across the pipe area and along the pipe length gives the rate of energy dissipation

Fig. 3. Control volume applied to the pipe for a standard test rig of water hammer.

129

A. Riasi et al. / Computers & Fluids 73 (2013) 124–133 Table 1 Grid independency study.

5. Results and discussion

Mesh grid

Time step (s)

385  47 500  62 600  80

0.00007 0.000053 0.000041

For every numerical works, two important issues should be checked: (a) numerical uncertainty (verification) and (b) comparison with experimental results (validation). These are fully described in below. 5.1. Verification

Mesh Grid:385 × 47, Time Step: 0.00007 S Mesh Grid:500 × 62, Time Step: 0.000053 S

1.5

Mesh Grid:650 × 80, Time Step: 0.000041 S 1

(H-H0)/(aV0 /g)

0.5 0 -0.5 -1 -1.5

0

2

4

6

8

at/L Fig. 4. Effect of grid size on pressure oscillations at the valve.

/0 ¼

Z

L

0

Z

R

/ 2prdr dx

ð25Þ

o

The energy dissipation in transient pipe flow may be computed by integrating Eq. (25) during a computational time interval, T, as

/00 ¼

Z

T

/0 dt

The main purpose of verification (including the code and solution verification) is to check the numerical accuracy independent of the physical modeling which is the subject of validation [15]. Regarding the solution verification of the numerical scheme, a standard test rig data of Holmboe and Rouleau [16] has been considered. Three different meshes including a coarse grid of 385  47 points, a fine mesh with 500  62 grids nodes and a finer mesh with 650  80 grid nodes, has been performed (Table 1). In all of these meshes, the clustering parameter is set to 1.05 in the radial direction. The main characteristics of the unsteady turbulent flow including pressure, turbulent velocity profile and kinetic energy of turbulence are considered. The comparison among the results obtained from these different meshes shown in Figs. 4–6. The details of uncertainty estimation based on the grid convergence index (GCI) are calculated and listed in Table 2 [15]. The numerical uncertainties are in the acceptable margin. where U1, U2 and U3 are values of key variables for finer, fine and coarse meshes, respectively; e21 = U2  U1, e32 = U3  U2, p is observed order of method, e21 is error in non-dimension form, GCI denotes the grid convergence index and unum is standard uncertainty in the simulation solution [15]. The following variables are considered for the key variables: (a) the maximum deviation for the valve pressure, (b) the maximum deviation for the steady state turbulent velocity and the steady state kinetic energy, (c) the maximum deviation for the unsteady turbulent velocity and the unsteady kinetic energy at the middle of pipe and at t = 4.5L/a.

ð26Þ

0

Using Eqs. (22) and (26), one obtains

5.2. Validation

Z Z /00 ¼ uqdVol  uqdVol c:v : c: v : t¼T t¼0 ! ! Z Z V2 V2 ¼ qdVol  qdVol 2 2 c:v : c:v :

To validate the unsteady friction model developed in this paper, the computed data are compared with the measured data of Holmboe and Rouleau [16] for the laminar and the turbulent water hammer. When the downstream valve is suddenly closed, a strong transient flow is produced in the pipe line (Fig. 3).

t¼0

ð27Þ

t¼T

Mesh Grid:385 × 47, Time Step: 0.00007 S Mesh Grid:500 × 62, Time Step: 0.000053 S

Mesh Grid:385 × 47, Time Step: 0.00007 S Mesh Grid:500 × 62, Time Step: 0.000053 S Mesh Grid:650 × 80, Time Step: 0.000041 S

1

Mesh Grid:650 × 80, Time Step: 0.000041 S

1 0.8

0.6

0.6

r/R

r/R

0.8

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

0

0

0.2

0.4

0.6

0.8

u/V0

u/V0

(a)

(b)

Fig. 5. Effect of grid size on u-velocity profile: (a) steady state and (b) t = 4.5L/a.

1

1.2

1.4

130

A. Riasi et al. / Computers & Fluids 73 (2013) 124–133 Mesh Grid:385 × 47, Time Step: 0.00007 S Mesh Grid:500 × 62, Time Step: 0.000053 S Mesh Grid:650 × 80, Time Step: 0.000041 S

4

4

3.5

3.5

3

3

2.5

2.5 2

k/uτ

k/uτ

2

2 1.5

2 1.5

1

1

0.5

0.5

0

0

-0.5

0.2

0

0.4

0.6

0.8

-0.5

1

Mesh Grid:385 × 47, Time Step: 0.00007 S Mesh Grid:500 × 62, Time Step: 0.000053 S Mesh Grid:650 × 80, Time Step: 0.000041 S

0

0.2

0.4

0.6

r/R

r/R

(a)

(b)

0.8

1

Fig. 6. Effect of grid size on kinetic energy of turbulence profile: (a) steady state and (b) t = 4.5L/a.

Table 2 The details of uncertainty estimation [15].

U1 U2 U3

e21 e32 p e21 GCI unum

Table 3 The detail of test cases examined.

Valve pressure

Steady state turbulent velocity

Unsteady turbulent velocity

Steady state kinetic energy

Unsteady kinetic energy

1.015 1.01226 1.00833 0.00274 0.00393 1.3747 0.27% 0.7769% 0.3884%

1.1795 1.1848 1.1908 0.0053 0.0059 0.4194 0.45% 4.8558% 2.4279%

0.7327 0.7443 0.7597 0.0116 0.0154 1.0636 1.58% 6.1744% 3.0872%

3.6466 3.5577 3.4338 0.0888 0.1239 1.2699 2.43% 7.6999% 3.8499%

2.0209 2.0627 2.1197 0.0418 0.0569 1.1789 2.06% 7.1355% 3.5677%

Test case number Reynolds Wave speed P (m/s) 1 2 3

6166 6166 6166

To study the influence of non-dimensional groups a and b on characteristics of turbulent water hammer problems, three different test cases with various values for a and b are examined. The details of these test cases are presented in Table 3. Note that the variations of a and b are achieved by changing the pressure wave velocity (we consider a very wide range, far behind

Exprimental (Holmboe et al., 1967) Computed

Computed Input Data for Laminar water Hammer: L=36.09 m D=0.025 m 3 ρ=878.4 Kg/m 2 μ=0.03484 N-s/m a=1324 m/s Re=82

1.5

0.5 0 -0.5 -1 -1.5 -2

Input Data for Turbulent water Hammer: L=36.09 m D=0.025 m 3 ρ=1000 Kg/m μ=0.00086 N-s/m2 a=1350 m/s Re=6166

2

(H-H0)/(aV0 /g)

(H-H0)/(aV0 /g)

1

4.75 4.33  103 0.269 47.5 4.37  104 0.0269 237.53 8.47  105 0.00538

135 1350 6750

5.3. Results discussion

Exprimental (Holmboe et al., 1967)

1.5

1 0.5 0 -0.5 -1 -1.5

0

2

4

6

8

10

12

b ¼ 0:01n DL M

are 0.0535 and 0.0742 for Fig. 7a and b, respectively. Since, the standard uncertainty in the experimental data is unknown; it is not possible to estimate the validation standard uncertainty. Regarding Fig. 7, the events occurs in a working cycle of water hammer (T = 4L/a) are fully described in the Section 3.

Fig. 7 presents the comparison of the computed and measured pressure head at the middle of the pipe and also for both a laminar and a turbulent water hammer. The initial data for these experiments is shown in this figure as well. The time axis is normalized by the wave travel time from the reservoir to the valve (L/a). The pressure axis is also normalized by Joukowsky head (aV 20 =g). The validation comparison error, E, which is defined as difference between simulation solution value and experimental data values,

2

M a ¼ n DL Re

14

-2

0

2

4

6

8

at/L

at/L

(a)

(b)

10

12

14

Fig. 7. Pressure–time history at the pipe midpoint: (a) laminar water hammer, Re = 82 and (b) turbulent water hammer, Re = 6166.

131

A. Riasi et al. / Computers & Fluids 73 (2013) 124–133

α=8.47×10-5, β=0.00538

-5

α=8.47×10 , β=0.00538

-4

α=4.37×10 , β=0.0269

30

-4

α=4.37×10 , β=0.0269

α=4.37×10-3, β=0.269

-3

α=4.37×10 , β=0.269

1

25 20

φ' /φ' ss

(H-H0)/(aV0 /g)

0.5

0

15 10

-0.5

5 0

-1 0

2

4

6

8

10

12

0

2

4

6

at/L Fig. 8. Pressure–time history at the pipe midpoint for various values of a and b.

-4

-3

α=4.37×10 , β=0.269

30 First step

Test case number

t (s)

U00v is (J)

U00tur (J)

1 2 3

3.2 0.32 0.064

0.285 0.1185 6.61  102

1.44  102 5.335  103 1.79  103

Second step

10

τW /τWss

12

Table 4 Energy dissipation in a three cycle of water hammer.

α=4.37×10 , β=0.0269

ratio of the diffusion time scale to the wave time scale in water hammer problems and defined as [17]

0

P¼

-10 -20 -30 -40

10

Fig. 11. Energy dissipation time history at the pipe midpoint for various values of a and b.

α=8.47×10-5, β=0.00538

20

8

at/L

0

2

4

6

8

10

12

at/L Fig. 9. Wall-shear-stress time history at the pipe midpoint for various values of a and b.

the limits that can be found in practice for only a parametric study). Case (2) represents the experimental data of Holmboe and Rouleau [16]. Also, the non-dimensional parameter P is the

pffiffiffi 2D=2us L=a

ð28Þ

Fig. 8 shows the head pressure at the middle of the pipe for the cases (1–3). It is evident that the pressure damping is much stronger for test case (1) that the values of a and b are larger than others. Two asymptotic limits of the momentum equation can be distinguished as: (a) When a, b ? 0 and as a result the viscous and turbulent terms are neglected and the pressure damping tends to be zero. On the other hand, the condition of inviscid water hammer is achieved. (b) When a, b ? 1 and consequently the unsteady and pressure terms are neglected and the period of the transient flow tends to infinity and as a result the transient nature of problem is vanished in the recent condition.

1

1

0.95

0.95

0.95

0.9

0.9

0.9

0.85

r/R

1

r/R

r/R

( ) Steady State, ( ) t=0.5L/a, ( ⋅ ) t=1.5L/a, ( ⋅⋅ ) t=2.5L/a

0.85

0.85

0.8

0.8

0.8

0.75

0.75

0.75

0.7

0.7

-1.5

-0.5

0.5

1.5

-1.5

-0.5

0.5

1.5

0.7

-1.5

-0.5

u/V0

u/V0

u/V0

(a)

(b)

(c)

0.5

1.5

Fig. 10. Unsteady velocity profiles for various values of a and b: (a) a = 4.37  103, b = 0.269, (b) a = 4.37  104, b = 0.0269 and (c) a = 8.47  105, b = 0.00538.

132

A. Riasi et al. / Computers & Fluids 73 (2013) 124–133

0.15

0.1

0.6

0.8

1

-0.1

-0.05

0

0.2

0.4

0.6

0.8

1

-0.1 -0.15

-0.15

0.1

(a)

0 4

0.4

0.15

0.05

r/R νε/u τ

0.2

4

0

0

r/R νε/u τ

4

νε/u τ

-0.05

0.1 0.05

0.05 0

0.2 4

0.15

νP/uτ

0.2 4

0.2

νP/uτ

νP/uτ

4

(−) Steady State, (--) t=0.5L/a, (-⋅-) t=1.5L/a, (-⋅⋅-) t=2.5L/a, (− −) t=3.5L/a

0.2

0.4

0.6

0.8

r/R

-0.05 -0.1 -0.15

(b)

(c)

Fig. 12. Profiles of turbulence production and turbulence dissipation for various values of a and b: (a) a = 4.37  103, b = 0.269, (b) a = 4.37  104, b = 0.0269 and (c) a = 8.47  105, b = 0.00538.

The oscillations of the wall shear stress are shown in Fig. 9. According to this figure, two different steps can be identified for the wall shear stress oscillation. The first step is a sudden rise of wall shear stress which is caused directly by the passage of the pressure wave. Kucienska [18] showed this shear stress peak is approximately proportional to the amplitude of the pressure wave. The second step is the relaxation of the wall-shear-stress. During this relaxation phase, the shear stress tends exponentially from the value caused by the pressure wave to the value corresponding to the new steady state. This situation is achieved if the flow conditions remained unchanged after the passage of the pressure wave which had caused the shear stress peak. The sudden pressure rises of wall shear stress are caused by sudden pressure gradients and they are responsible for accelerating and decelerating flows. Fig. 10 shows the velocity profiles at the first wave cycle of water hammer. It is clearly shown when the positive pressure wave passes from the pipe midpoint (t = 0.5L/a), a uniform shift in the velocity profile is induced in the whole of pipe cross-section. In addition, because of the no-slip boundary condition at the pipe wall, a strong gradient and a vortex ring (reverse flow) is formed around the pipe. The same condition occurs when the negative pressure wave coming from the valve (t = 2.5L/a). During deceleration of flow, the velocity and the shear stress near the wall have signs opposite to the velocity in the core flow and reverse flow occurs near wall. This phenomenon called annular effect and was first introduced by Richardson and Tyler in 1929. Regarding Figs. 9 and 10, the sudden rise of wall shear stress and also annular effect for test case (3) is much significant when compared to those happens in test cases (1) and (2). Fig. 11 represents the rate of energy dissipation in three cycles of water hammer (12L/a). As shown in this figure, the maximum rate of energy dissipation occurs at the time t = 2L/a when the flow in the whole of the pipe is accelerated because of traveling of reflected positive pressure wave from the reservoir to the valve. Referring to this figure, the sudden rise of energy dissipation rate reduces as the a and b are increased. The turbulent and viscous energy dissipation (U00v is ; U00tur ), during the three wave cycles and for the mentioned test cases are listed in Table 4. Note that because of various wave speeds, the necessary time for three wave cycles is different in the three mentioned test cases. Fig. 12 indicates the profiles of turbulent production and dissipation in the first wave cycle of water hammer and at the middle of the pipe. The peaks of turbulent production and dissipation occur in the wall region where the velocity gradient is substantial and tend to zero at the pipe axis. For test case (1) that P  1 and the values of a and b are larger, the turbulent production and dissipation

are considerably changed from their steady state values and the peaks move away from the wall. This variations are small for case of P  1, for where the values of a and b are smaller. 6. Conclusions In this paper, the energy dissipation and turbulent production and dissipation in unsteady turbulent pipe flows due to water hammer problems, has been numerically investigated. For the both steady and unsteady turbulent pipe flow, a good agreement between the numerical and the experimental results is observed. The results show the k–x turbulence model has reasonable capability in prediction of unsteady pipe flows due to water hammer problems. Through dimensional analysis and order of magnitude, two non-dimensional groups are identified that can be used for modeling the viscous and turbulent shear stresses. The result shows that when a and b decrease: (a) the pressure damping decreases, (b) the sudden rise of wall shear stress and dissipation function increase and (c) the effect of Richardson annular flow is more pronounced. Moreover, the non-dimensional parameter, P, decreases as the values of a and b increase. In case of P  1 and for a cycle of water hammer, the turbulence structure is significantly changed and the peaks of turbulent production and dissipation move rapidly from the wall to the pipe axis. This variations are small for case of P  1, for which the values of a and b are smaller. Acknowledgment The authors wish to extend their utmost gratitude to the staff of Hydraulic Machinery Research Institute (HMRI) for their cooperation and support. References [1] Silva-Araya WF, Chaudhry MH. Computation of energy dissipation in transient flow. J Hydraul Eng 1997;123(2):108–15. [2] Wahba EM. Modelling the attenuation of laminar fluid transients in piping systems. Appl Math Model 2008;32:2863–71. [3] Vardy AE, Hwang KL. A characteristic model of transient friction in pipes. J. Hydraul Res 1991;29(5):669–85. [4] Silva-Araya WF, Chaudhry MH. Unsteady friction in rough pipes. J Hydraul Eng 2001;127(7):607–18. [5] Pezzinga G. Evaluation of unsteady flow resistances by quasi-2D or 1D models. J Hydraul Eng 2000;126(10):778–85. [6] Zhao M, Ghidaoui MS. Efficient quasi two dimensional model for water hammer problems. J Hydraul Eng 2001;129(12):1007–13. [7] Wahba EM. Runge–Kutta time-stepping schemes with TVD central differencing for the water hammer equations. Int J Numer Methods Fluids 2006;52(5):571–90. [8] Zhao M, Ghidaoui MS. Investigation of turbulent behavior in pipe transient using a k–e model. J Hydraul Res 2006;44(5):682–92.

A. Riasi et al. / Computers & Fluids 73 (2013) 124–133 [9] Fan S, Lakshminarayana B, Barnett M. Low-Reynolds number k–e model for unsteady turbulent boundary layer flows. AIAA J 1993;31(10):1777–84. [10] Riasi A, Nourbakhsh A, Raisee M. Unsteady turbulent pipe flow due to water hammer using k–x turbulence model. J Hydraul Res 2009;47(4):429–37. [11] Tannehill JC, Anderson DA, Pletcher RH. Computational fluid mechanics and heat transfer. USA: Taylor and Francis; 1984. [12] Wilcox DC. Turbulence modelling for CFD. La Canada CA: DCW Industries; 1994. [13] Hinze JO. Turbulence. 2nd ed. New York: McGraw-Hill; 1987. [14] White FM. Viscous fluid flow. 2nd ed. New York: McGraw-Hill; 1991.

133

[15] Standard for verification and validation in computational fluid dynamics and heat transfer. The American Society of Mechanical Engineering; 2009. [16] Holmboe EL, Rouleau WT. The effect of viscous shear on transients in liquid lines. Trans ASME 1967;89(1):174–80. [17] Ghidaoui MS, Mansour SGS, Zhao M. Applicability of quasisteady and axisymmetric turbulence models in water hammer. J Hydraul Eng 2002;128(10):917–24. [18] Kucienska B. Friction relaxation model for fast transient flows. PhD thesis. Belgume: University of Catholique de Louvain; 2004.