On the dissipation mechanism of lattice Boltzmann method when modeling 1-d and 2-d water hammer flows

Accepted Manuscript

On the dissipation mechanism of lattice Boltzmann Method when modeling 1-d and 2-d water hammer flows Moez Louati, Mohamed Mahdi Tekitek, Mohamed Salah Ghidaoui PII: DOI: Reference:

S0045-7930(18)30584-X https://doi.org/10.1016/j.compfluid.2018.09.001 CAF 3996

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

5 April 2017 25 July 2018 1 September 2018

Please cite this article as: Moez Louati, Mohamed Mahdi Tekitek, Mohamed Salah Ghidaoui, On the dissipation mechanism of lattice Boltzmann Method when modeling 1-d and 2-d water hammer flows, Computers and Fluids (2018), doi: https://doi.org/10.1016/j.compfluid.2018.09.001

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • All requests of reviewers are made.

CR IP T

• The writing is improved with the help of native English expert in transient flow.

• The remarks and advice of review made the paper clearer and more con-

AC

CE

PT

ED

M

AN US

sistent.

1

ACCEPTED MANUSCRIPT

Moez Louati1

CR IP T

On the dissipation mechanism of lattice Boltzmann Method when modeling 1-d and 2-d water hammer flowsI

Department of Civil and Environmental Engineering, The Hong Kong University of Science and Technology.

Mohamed Mahdi Tekitek

Department of Mathematics, Faculty of Science of Tunis, University of Tunis El Manar, 2092, Tunis, Tunisia.

AN US

Mohamed Salah Ghidaoui

Department of Civil and Environmental Engineering, The Hong Kong University of Science and Technology.

M

Abstract

This paper uses the multiple relaxation times Lattice Boltzmann Method (LBM)

ED

to model acoustic wave propagation in water-filled conduits (i.e. water hammer applications). The LBM scheme is validated by solving some classical water hammer (WH) numerical tests in one-dimensional and two-dimensional cases.

PT

In addition, this work discusses the accuracy, stability and robustness of LBM when solving high frequency (i.e. radial modes are excited) acoustic waves in water-filled pipes. The results are compared with a high order finite volume

CE

scheme based on Riemann Solver. The results show that LBM is capable to model WH applications with high accuracy and performance at low frequency

AC

cases; however, it loses stability and performance for high frequency (HF) cases. It is discussed that this is due to the neglected high order terms of the equivalent

I Contribution submitted for publication, Proceedings of the ICMMES Conference, July, 2016. Email addresses: [email protected] (Moez Louati ), [email protected] (Mohamed Mahdi Tekitek), [email protected] (Mohamed Salah Ghidaoui) 1 corresponding author

Preprint submitted to Computers and Fluids

July 24, 2018

ACCEPTED MANUSCRIPT

macroscopic equations considered for LB numerical formulation. Keywords: LB scheme; Water hammer; MRT; Numerical dissipation; High

CR IP T

frequency; Acoustic waves.

Introduction

Waves are used as diagnostic tools for imaging and defect detection in a

wide range of science and engineering applications (e.g. medicine, oceanography, materials and many others. In particular, acoustic waves are used in sonog-

AN US

raphy (echography in medicine) and ocean topography identification. Current work in hydraulic (engineering) research uses acoustic waves for defect detec-

tion in water supply systems (e.g. [13, 29, 34, 36]). More recently, attempts are being made to spatially locate and identify defects at high resolution using high frequency waves (HFW). However, the use of high frequency probing waves (HFW) excites radial and azimuthal wave modes in pipe systems, making

hydraulic problem.

M

classical one-dimensional water hammer (WH) theory incapable of solving the

To conserve the physical dispersion behavior of HFW, dispersive numeri-

ED

cal schemes (e.g. [31, 32, 33]) should be avoided. Louati and Ghidaoui [33] showed that a 2nd order FV scheme performs poorly for modeling HFW and

PT

high order schemes are required to reduce numerical dispersion effects. Louati and Ghidaoui [31, 32] used a fifth order numerical scheme based on a Riemann solver with WENO reconstruction to model and study the behavior of HFW

CE

in a simple pipe system. However, because of complex boundary conditions needed in hydraulic engineering (e.g., pumps, pressure relief valves, pressure

AC

reducing valves, and so on), it is best to avoid the use of high order schemes due to the increased mathematical and numerical difficulty of implementing the boundary conditions. The lattice Boltzmann method (LBM) is a popular alternative to high order schemes and is known for its excellent performance in engineering applications. For this reason, this paper studies the efficiency of the LBM to model water hammer flows and high frequency waves. In fact, LBM

2

ACCEPTED MANUSCRIPT

makes it possible to simulate various types of fluid flows with simple algorithms (see e.g. [7, 11, 28, 42, 43]). In simple cases, second order accuracy can be easily verified (see e.g. [27]). As the LBM needs only nearest neighbor infor-

CR IP T

mation to achieve a solution, the algorithm is an ideal candidate for parallel computing. These features make the LBM increasingly popular for engineering applications. Currently, the LBM is a well-known computational fluid dynamics method that models flow in a way that is consistent with formal solutions of

viscous and incompressible Navier-Stokes equations [40]. However, use of the LBM for acoustics (compressible flow) applications is fairly recent and its utility

AN US

and performance characteristics are not yet well known.

LBM was first used by Buick et al. [4] to simulate sound waves (at low frequency) in situations where fluid density variation is small compared with mean density. Dellar [10] studied the ability of the LB scheme to model sound waves and investigated the influence on the numerical scheme of bulk viscosity and shear viscosity. Course et al. [9] used LBM to investigate four canonical

M

acoustic problems involving traveling and standing acoustic waves and showed the capability of the LBM to reproduce fundamental acoustic wave phenomena.

ED

Recently Mari´e et al. [35] showed that dispersion error of acoustic waves in lattice Boltzmann models is due to effects of space and time discretization. They also studied dispersion and dissipation errors of the LB scheme and compared

PT

the result to various Navier-Stokes solvers. Their results highlighted the low dissipative nature of lattice Boltzmann models compared with other high order

CE

schemes.

More recently, Budinski [3] investigated solving the 1-d water hammer equa-

tion, with emphasis on models that configure the computational grid indepen-

AC

dently of the wave speed. He proposed to use LBM method for one-dimensional problems using the D1Q3 model with three discrete velocities. The governing equations are modified using appropriate coordinate transformations to eliminate grid limitations related to the method of characteristics. This paper addresses the following key issues :

3

ACCEPTED MANUSCRIPT

1. Modeling water hammer: In the field of hydraulic engineering and transient flows, LBM is not considered to be suitable for modeling compressible flow, and generally the method is avoided. The method of characteristics

CR IP T

(MOC) is preferred for modeling transient flows in 1-d and finite volume method or finite difference method are typically used for 2-d water ham-

mer flows. However, this paper shows that LBM is equally accurate as

MOC and, since LBM is computationally a very fast method, it could be used to efficiently model 1-d transient flows. Therefore, in the first part of this paper, classical 1-d water hammer applications (e.g. Ghidaoui [22],

AN US

Chaudhry [6]) are modeled using LBM and compared with MOC results to

evaluate the accuracy and performance of LBM scheme relative to MOC. 2. In the second part of this paper, the 2-d water hammer case is modeled using LBM. Comparison with a 5th order finite volume (FV) scheme shows that LBM suffers from dissipation and dispersion at sharp discontinuities and when simulating radial oscillations. Notice that the preservation of

M

sharp and continuous wavefronts in water hammer actually requires maintenance of longitudinal waves at high frequency in the numerical scheme.

ED

Although a sharper wavefront could be maintained in LBM by mesh refinement, an alternative approach is considered in this paper by investigating the dissipation effects in LBM considering smooth high frequency waves

PT

where the longitudinal (plane wave) mode and radial modes propagate separately.

CE

Because radial and longitudinal wave modes may have different dissipation

mechanisms, multi-relaxation time (MRT) LBM is used in this study to allow separate treatment of the bulk and shear viscosity terms. In fact, the results

AC

show that radial wave dissipation is governed by the bulk viscosity whereas dissipation of longitudinal waves is dictated by shear viscosity. In addition, the MRT scheme has additional free parameters that have no analogous physical interpretation but which do affect the stability of the LBM scheme and also affect accuracy of the boundary conditions [21] solutions. Other variants of

4

ACCEPTED MANUSCRIPT

the LBM scheme can also be used, such as two relaxation times (TRT) [24], single relaxation time (SRT) [24], or other LBM schemes using different collision operators such as entropic lattice Boltzmann (ELBM) [8] and cumulant lattice

CR IP T

Boltzmann (CLBM) [23]. In fact, the ELBM scheme claims that it respects the discrete H theorem while CLBM uses more general equilibrium distributions

than those of the MRT scheme. The choice to use MRT LBM in this work is

based on its simplicity and popularity. Future work will examine in more depth the performance of other LBM schemes.

AN US

1. Lattice Boltzmann scheme for water hammer application 1.1. One-dimensional case

This section describes the lattice Boltzmann scheme for one-dimensional water hammer application. The LB scheme uses three discrete velocities (D1Q3) as shown in Figure 1. Let L be a regular lattice parametrized by a space step

M

∆x (see Figure 1), composed by a set L0 ≡ {xj ∈ (∆xZ)} of nodes (or vertices); ∆t is the time step of the evolution of LB scheme and λ ≡

∆x ∆t

is the numerical

ED

celerity. The discrete velocities vi , i ∈ (0, 1, 2) are chosen such that vi ≡ ci ∆x ∆t =

ci λ, where the family of vectors {ci } is defined by: c0 = 0, c1 = 1 and c2 = −1. The lattice Boltzmann scheme is given as follows [14]:

PT

fi (xj , t + ∆t) = fi∗ (xj − vi ∆t, t),

0 ≤ i ≤ 2,

(1)

where fi is the discrete probability distribution function corresponding to the

CE

ith discrete velocity vi , t denotes time, xj is the position of the j th node and the superscript ∗ denotes post-collision quantities. Therefore, during each time

AC

increment ∆t there are two fundamental steps: advection and collision. The advection step describes the motion of the particles traveling from nodes xj −∆x and xj + ∆x with velocities v1 and v2 , respectively, towards node xj where the

collision step occurs [25]. The collision process is defined by the conservation of moments of order zero (mass) and of order one (momentum) from the precollision state to the post-collision state between advected particles and the local 5

ACCEPTED MANUSCRIPT

particle at node xj [25]. The three moments {m` , ` ∈ (0, 1, 2)} are obtained by a linear transformation from the space of discrete probability density functions. These moments are density m0 ≡ ρ ≡ f0 + f1 + f2 , momentum m1 ≡ q ≡ ρu ≡ λ2 2 (f1

+ f2 ). The

CR IP T

λ(f1 − f2 ) (where u is the fluid velocity) and energy m2 ≡  ≡

linear transformation is summarized in matrix form as follows: (m0 , m1 , m2 )t = M (f0 , f1 , f2 )t , where



1

0

λ 2

λ 2

0

1



  −λ  .  2

λ 2

AN US

  M = 

1

(2)

(3)

The non-conserved moment (energy m2 ) is assumed to relax towards the equilibrium value (meq 2 ) as follows [15]:

m∗2 = (1 − s) m2 + s meq 2 ,

(4)

where s (s > 0) is the relaxation rate and the equilibrium value meq 2 is defined

M

as


1 αλ2 ρ + ρ u2 . (5) 2 R 0 2 The above expression Eq. (5) is obtained from the second order moment R f u du

ED

meq 2 =

where f 0 is the Maxwellian distribution function. Consequently, using Taylor

expansion at order 2 [15], the following two macroscopic equations are obtained

AC

CE

PT

from the advection and collision steps :   ∂t ρ + ∂x q = O(∆t2 ),  ∂ q + c2 ∂ ρ + ∂ (ρu2 ) − ν ∂ 2 q = O(∆t2 ), t x 0 x s x √ where cs = λ α is the equivalent sound speed and   1 1 − ν0 = λ2 ∆t(1 − α) s 2

(6)

(7)

is the equivalent kinematic viscosity. The meq 2 moment Eq. (5) determines the macroscopic behavior Eq. (6) of the scheme Eq. (1). The state equation relating the pressure to the fluid density is ∂P = c2s , ∂ρ 6

(8)

ACCEPTED MANUSCRIPT

where cs is the wave speed. Eq. (8) is for rigid pipe and in this case cs is about 1440 m/s when the fluid is water. However, for simplicity in this work, cs is taken equal to 1000 m/s although Eq. (8) is used.

CR IP T

Note that in this one-dimensional model, there is only one relaxation parameter (s). Therefore, with the linear transformation between the {m` } and {fi } spaces, Eq. (2) and the equilibrium distribution in moment space : m0 = meq 0 = ρ,

m1 = meq 1 = q,

meq 2 =


1 αλ2 ρ + ρu2 , 2

the following thermodynamic equilibrium is obtained (in f space):

=

eq 2 λ2 m 2

ρ u2 , λ2

AN US

2

= meq 0 −

1 eq 1 m + 2 meq 2λ 1 λ 2

1 1 eq = − meq 1 + 2 m2 2λ λ

=

(1 − α)ρ −

=

ρu 1 + 2λ 2

  ρ u2 αρ + 2 , λ

ρu 1 = − + 2λ 2

M

   f0eq           f1eq             f eq



ρ u2 αρ + 2 λ



(9)

.

ED

which describes the classical BGK distribution (see [28]).

1.2. Two-dimensional case

PT

In the two-dimensional case, the regular lattice L becomes as shown in Fig-

ure 2 comprising a set L0 ≡ {xi,j ∈ (∆xZ, ∆yZ)} of nodes (or vertices) with ∆x = ∆y. Similar to the previous one-dimensional case, the LB scheme in

fi (xj , t + ∆t) = fi∗ (xj − vi ∆t, t),

0 ≤ i ≤ 8,

(10)

AC

CE

two-dimensions (2-d) is given by:

which in this case, uses nine discrete velocities vi = λ ci (D2Q9) (see Figure

2) where c0 = (0, 0), c1 = (1, 0), c2 = (0, 1), c3 = (−1, 0), c4 = (0, −1), c5 = (1, 1), c6 = (−1, 1), c7 = (−1, −1) and c8 = (1, −1). Again, the scheme is based on the advection and collision steps as described in the 1-d case. However, in the 2-d case the scheme deals with nine moments {mk , k = 0, . . . , 8}. These 7

ACCEPTED MANUSCRIPT

moments have an explicit physical meaning (see e.g. [28]): m0 = ρ is the density, m1 = jx and m2 = jy are x−momentum, y−momentum, m3 is the energy, m4 is proportional to energy squared, m5 and m6 are x−energy and

CR IP T

y−energy fluxes and m7 , m8 are diagonal and off-diagonal stresses [28]. The same linear transformation Eq. (2) is used to obtain the domain of moments from the domain of distribution functions {fk , k = 0, . . . , 8} where the matrix M is given in Eq. (A.1). The equivalent Navier-Stokes equations are obtained by conserving three moments namely: m0 = ρ density, m1 = jx = ρVx and

m2 = jy = ρVy the x−momentum and y−momentum, respectively. The other

AN US

non-conserved moments are assumed to relax towards equilibrium values (meq ` ) obeying the following relaxation equation [25] : m∗` = (1 − s` ) m` + s` meq ` ,

3 ≤ ` ≤ 8,

(11)

where the equilibrium values are given by Eq. (A.2) and where s` (0 < s` < 2, for ` ∈ {3, 4 . . . , 8}) are relaxation rates, not necessarily equal to a single

M

value as in the so-called BGK scheme [40]. The coefficients α and β are tuning parameters which will be fixed later.

ED

Expanding Eq. (10) using Taylor series, the equivalent macroscopic equations of order two are Eq. (A.3), (A.4) and (A.5).

PT

Inserting jx = ρVx and jy = ρVy in the equivalent macroscopic equations (see

∂ρVy ∂ρ ∂ρVx + + ∂t ∂x ∂y

AC

CE

Eq. (A.3), Eq. (A.4) and Eq. (A.5)), yields:

8

=

0,

(12)

ACCEPTED MANUSCRIPT

     λ2 −α 1 1 ∂ ∂Vx ∂Vy ∆t − ρ + + 3 2 s3 2 ∂x ∂x  ∂y     2 1 ∂2 1 ∂ + + 2 Vx − ρ s8 2 ∂x2 ∂y

+



+

λ2 ∆t 3

+

λ2

=



CE +



1 1 − s3 2

1 1 − s8 2

 ∂ρ ∂Vx ∂2ρ Vx 2 + 2 + ∂x ∂x ∂x  2 ∂ρ ∂Vy ∂ρ ∂Vy ∂ ρ + + +Vy ∂x∂y ∂x ∂y ∂y ∂x

  ∂2ρ ∂ρ ∂Vx ∂2ρ ∂ρ ∂Vx + Vx 2 + 2 Vx 2 + 2 , ∂x ∂x ∂x ∂y ∂y ∂y (13)

α + 4 ∂ρ ∂ ∂ + ρVx Vy + ρV 2 = 6 ∂y ∂x ∂y y

     λ2 −α 1 ∂ ∂Vx ∂Vy 1 ∆t − ρ + 3 2 s3 2 ∂y ∂x  ∂y   2   1 ∂ 1 ∂2 + − ρ + 2 Vy s8 2 ∂x2 ∂y 

−αλ2 ∆t 6

PT

+

−αλ2 ∆t 6

M

∂ρVy ∂t

α + 4 ∂ρ ∂ ∂ + ρV 2 + ρVx Vy = 6 ∂x ∂x x ∂y

CR IP T

=

λ2

AN US

+

ED

∂ρVx ∂t

λ2 ∆t 3





1 1 − s3 2

1 1 − s8 2

 ∂ρ ∂Vy ∂2ρ + Vy 2 + 2 ∂y ∂y ∂y  ∂ρ ∂Vx ∂ρ ∂Vx ∂2ρ + + Vx ∂x∂y ∂x ∂y ∂y ∂x

  ∂2ρ ∂ρ ∂Vy ∂2ρ ∂ρ ∂Vy Vy 2 + 2 + Vy 2 + 2 ∂x ∂x ∂x ∂y ∂y ∂y

(14)

AC

where Vx and Vy are the flow velocities along the x and y directions, respectively. Notice that the above equivalent macroscopic equations contain additional terms (shown between curly brackets in Eq. (13) and Eq. (14)) with respect to the Navier-Stokes equations. It can be shown that these additional terms are negligible with respect to the viscous terms.

9

ACCEPTED MANUSCRIPT

1.3. 2-d classical water hammer test case This section considers the case of viscous-laminar flow as shown in Figure 3. Notice that in this case y represents the radial coordinate and Vy is the radial

CR IP T

velocity. The initial conditions are constant initial velocity and pressure along the pipe given as follows:    Vx (y)         P (x)

=

2Vx0

= p0 −



y2 1− 2 R

32 ρ0 Vx0 D2



ν

,

x,

0 ≤ y ≤ R,

(15)

0 ≤ x ≤ L.

AN US

where ν = kinematic viscosity in water (ν = 10−6 m2 /s); D = 2R = pipe

diameter (D = 0.4 m) with R the pipe radius; Vx0 = constant mean velocity; ρ0 = initial density ( ρ0 = 1000 kg/m3 ); p0 = constant initial gauge pressure in the pipe. It is important to note that, in a cylindrical coordinate system, additional geometrical terms are added to the 2-d Navier-Stokes equations which

M

are given below:

∂ρVx ∂ρ ∂ρVy + + ∂t ∂y ∂x

ρVy = − y | {z }

(16)

ED

Term due toradial coordinates system

CE

PT

∂ρVy ∂ ∂ ∂ρ + c2s + ρVy2 + ρVy Vx ∂t ∂y ∂r ∂x

AC

∂ρ ∂ ∂ ∂ρVx + c2s + ρVy Vx + ρV 2 ∂t ∂x ∂y ∂x x

    ∂ ∂Vy ∂Vx = (µ + λ) + + µ∆Vy ∂y ∂y ∂x +

2µ + λ y |



∂Vy Vy − ∂y y {z



ρ − Vy2 y }

,

Terms due to radial coordinates system

(17)     ∂ ∂Vy ∂Vx = (µ + λ) + + µ∆Vx ∂x ∂y ∂x +

µ + λ ∂Vx µ ∂Vx ρ + − Vy Vx y ∂x y ∂y y {z } |

,

Terms due to radial coordinates system

10

(18)

ACCEPTED MANUSCRIPT

where Vx is the x−velocity and Vy is the y−velocity. The D2Q9 scheme introduced in the above section is used to solve Eq. (16), Eq. (17) and Eq. (18). Note that by using LBM, the equivalent equations are slightly different from the

CR IP T

continuous problem. The additional geometrical terms are added to the LBM scheme using operator splitting of order 2. 1.4. Boundary conditions

The boundary conditions are essentially to impose a non slip condition on the

wall or to impose a given pressure/velocity in upstream/downstream boundaries.

AN US

To perform these boundary conditions, anti-bounce back [24] scheme and bounce back [21] scheme are used to impose a given pressure and a given velocity, respectively. The detail of bounce back and anti-bounce back schemes are given in next section.

2. Numerical results and discussion

M

2.1. 1-d classical water hammer test case

The 1-d water hammer test case consists of sudden closure of a downstream

ED

valve in a reservoir-pipe-valve (RPV) system with an initial non-zero flow [44] (see Figure 4). The boundary conditions require pressure to be constant at the upstream (reservoir) for all time, and the velocity to be zero at the downstream

PT

(Vx = 0). These boundary conditions are imposed as follows: • Upstream boundary P (x = 0) = p0 : To impose pressure on the inlet, anti-

CE

bounce back [5] boundary condition is performed:

p0 c2s

(19)

is given by the boundary condition and x1 is the first fluid node.

AC

where ρ0 =

f1 (x1 ) = −f2∗ (x1 ) + αρ0 ,

Remark The above anti-bounce back condition is obtained by taking the distribution f1 and f2 at equilibrium (see Eq. (9)) at the boundary.

• Downstream boundary Vx (x = L) = 0: To model Vx = u0 = 0 on the outlet

bounce back boundary condition is performed: f2 (xN ) = f1∗ (xN ) + 11

ρu0 , λ

(20)

ACCEPTED MANUSCRIPT

where u0 = 0 is given by boundary condition and xN is the last fluid node. Remark To keep the stability of the scheme, the parameter α must be in the interval [0, 1] because the term (1 − α) in Eq. (6) must be positive.

which is defined by: α = cs

∆t . ∆x

AN US

√

CR IP T

For all cases, the length of the pipe is L = 1000 m and the mesh size is Nx = 101. √ Thus, the space step is ∆x = L/Nx and the numerical wave speed is cs = ∆x ∆t α. √ ∆x The time step and relaxation parameter are given by: ∆t = α and s = cs  −1 ν0 1 , respectively. Only the parameter α remains unfixed. + λ∆x(1 − α) 2 In fact this parameter represents the Courant–Friedrichs–Lewy (CFL) number,

(21)

Note that for the case α = 1 the equivalent kinematic viscosity (ν0 ) (see Eq. (7)) becomes zero, and therefore the LB macroscopic equations (see Eq. (6)) are equivalent to the inviscid 1-d water hammer equations [22].

Three test cases are considered for different values of α and s :

M

1. α = 1, in this case the numerical viscosity (ν0 ) is zero for any s value. Therefore the flow is inviscid.

ED

2. α = 1/3 and s is such that ν0 = 10−6 m2 /s which is equal to the actual kinematic viscosity value for water.

PT

3. α = 1/3 and s = 1.1.

Figures 5, 6 and 7 give dimensionless pressure variations with time measured at the valve. The pressure is normalized by the Joukowsky pressure (PJou =

CE

ρ0 Vx0 cs ) for cases 1, 2 and 3, respectively. Case 1 (α = 1) is depicted in Figure 5, which shows that LBM gives the exact solution of water hammer pressure

AC

variation for a RPV system with inviscid flow. The exact solution is known from the method of characteristics (MOC) solution and can be found in [44]. This result is interesting because LBM usually solves viscous flows. In addition the LB scheme (in the case α = 1 and s = 1) is equivalent to the

12

ACCEPTED MANUSCRIPT

following finite difference scheme:   ρn+1 = 1 ρn + ρn
+ 1 q n − q n
, j+1 j−1 j+1 j−1 j 2 2λ  q n+1 = 1 q n + q n
− λ ρn − ρn
, 2

j+1

j−1

2

j+1

(22)

j−1

CR IP T

j

which is the solution of the density and momentum obtained by the method of

characteristics; where n denotes the nth time step and j is the position of the j th node.

In case 2, viscous flow is considered, and viscosity is chosen to be equal to

the actual kinematic viscosity value for water (ν0 = 10−6 m2 /s). In this case,

AN US

Figure 6 shows that the scheme reaches its limit of stability. This is because

the parameter s governs the viscous term in Eq. (6) and becomes almost 2 (s = 1.999999999650126) whereas s must be strictly less than 2 for the scheme to remain stable.

In case 3, typical values of α and s are taken (α = 1/3 and s = 1.1). For these values, the numerical kinematic viscosity is about ν0 = 5000 m2 /s. However,

M

Figure 7 shows that the signal is slightly damped. Therefore, the numerical viscosity does accurately represent the actual flow viscosity. This may be due

Eq. (6)).

ED

to higher order terms neglected when deriving the macroscopic equations (see

PT

2.2. 2-d classical water hammer test case This section considers the problem described in the section 1.3, and solves the macroscopic equations Eq. (12), Eq. (13) and Eq. (14) using LB. Con-

CE

sider the domain Ω = [0, L] × [0, H] and the regular mesh (lattice) L0 = {xi,j ∈ ∆xZ×∆yZ, 1 ≤ i ≤ Nx , 1 ≤ j ≤ Ny } parametrized by the space step ∆x = ∆y.

AC

The boundary conditions are like those in the previous section with the additional constraints of zero radial velocity and no slip condition at the pipe walls. In more detail, the boundary conditions are as follows: • Upstream boundary (P (x = 0) = p0 ) : To impose p0 , anti-bounce back is performed as follows:

13

ACCEPTED MANUSCRIPT

 p0   f (x ) = −f3 (x(0,j) ) + 4−α−2β  18 c2s ,  1 (1,j) 4+2α+β p0 f5 (x(1,j) ) = −f7 (x(0,j−1) ) + 18 c2 ,  s   p0  f (x −f6 (x(0,j+1) ) + 4+2α+β 8 (1,j) ) = 18 c2

(23)

CR IP T

s

where p0 is given by the boundary condition and 1 ≤ j ≤ Ny .

• Downstream boundary (Vx (x = L) = 0): To impose Vx = 0, bounce back is performed as follows:     

f3 (x(Nx ,j) ) = f1 (x(Nx +1,j) ),

where 1 ≤ j ≤ Ny .

(24)

AN US

f6 (x(Nx ,j) ) = f8 (x(Nx +1,j+1) ),     f (x 7 (Nx ,j) ) = f5 (x(Nx +1,j−1) )

• Pipe wall boundary (Vy (y = R) = 0 ): To impose Vy = 0, bounce back is

M

performed as follow: Bounce back for xi,1 and xi,Ny where 1 ≤ i ≤ Nx :    f (x ) = f2 (x(i,0) ),   4 (i,1) f5 (x(i,1) ) = f7 (x(i−1,0) ),     f (x ) = f (x ). 6

ED

   f (x )   4 (i,Ny ) f7 (x(i,Ny ) )     f (x )

PT

8

8

(i,1)

(i,Ny )

(25)

(i+1,0)

=

f2 (x(i,Ny +1) ),

=

f5 (x(i−1,Ny +1) ),

=

f6 (x(i+1,Ny +1) ).

(26)

The following length scales L = 1000 m, D = 0.4 m and mesh size Ny = 20, 160 and 500, Nx = cs λ,

Ny L Dqare

∆t considered. The CFL number is given by C0 = cs ∆x =

CE

where cs = λ α+4 = 1000 m/s is the wave speed. Therefore the CFL q 6 α+4 number is C0 = 6 . Thus, the LB parameters are fixed as follows: α =

AC

6C02 − 4, ∆t =

∆x cs C0 ,

α = −2, and β = 1. The time step ∆t and space step ∆x

are fixed to have wave speed cs = 1000 m/s. The other LB relaxation (stability)

parameters are fixed by identification between equations Eq. (13), Eq. (14) and Eq. (17), Eq. (18) as follows:   λ2 1 1 µ0 = ∆t − ρ = ν0 × ρ 3 s8 2 14

⇒

s8 =

1 , + 21

3ν0 λ∆x

(27)

ACCEPTED MANUSCRIPT

µ0 −αλ2 = ∆t 3 6



1 1 − s3 2



ρ

⇒

s3 =

1 2

−

2 3α

1 

1 s8

−

1 2

.

(28)

Eq. (27) satisfies the kinematic viscosity, and Eq. (28) satisfies Stokes relation.

CR IP T

For water hammer problem, where ν0 = 10−6 m2 /s, the stability parameters are s3 = 199999953 and s8 = 1.99999986. In this case the LB scheme is unstable because s3 and s8 approach the limit of stability (s3 < 2, s8 < 2). In what

follows, different tests are conducted to find the appropriate relaxation (stability) parameters for the 2-d water hammer application. Figures 8, 9, 10 and 11 give the dimensionless pressure variation with time measured at the valve

AN US

and at the pipe centerline where the pressure is normalized by the Joukowsky pressure (PJou = ρ0 Vx0 cs ) for the different cases. Figure 8 shows the case for the maximum s8 value to give a stable scheme while the Stokes relation Eq. (28) is maintained. In this case, the scheme is extremely dissipative, and maintain-

ing the Stokes relation is no longer possible. Figure 9 depicts the case where s8 = 1.99965 which is chosen such that ν0 = 10−3 m2 /s and s3 is given a typical

M

value of 1.1. In this case the signal becomes much less dissipative and gives satisfactory results. Note that the oscillations observed in Figures 8 and 9 are

ED

physical, and are due to radial waves (see [30, 31, 32]). These radial waves are excited by the idealized (instantaneous) valve closure and propagate along the pipe. However, Figure 9 shows that these waves are dissipated quickly. Figure

PT

10 gives the case where s8 is kept the same but s3 is increased to its maximum value (s3 = 1.45) while maintaining a stable scheme. In this case, Figure 10

CE

shows that the radial waves undergo less dissipation while the axial wave dissipation is unchanged. Figure 11 depicts the case where s3 is kept to its maximum value 1.45 while s8 is reduced to s8 = 1.9965 (ν0 = 10−2 m2 /s). In this case,

AC

Figure 11 shows that radial wave dissipation remains unchanged while the axial waves become more dissipated than in the previous two cases. This shows that s3 , which relates to the bulk viscosity, governs the dissipation mechanism of radial waves; whereas s8 governs the dissipation mechanism of axial waves. The values of s3 = 1.45 and s8 = 1.99965 are used for the rest of this paper.

Before continuing further discussion and comparison of the LB scheme with 15

ACCEPTED MANUSCRIPT

other numerical approaches, it is incumbent to conduct a convergence test using mesh refinement with the fixed relaxation parameters (s3 , s8 ). In fact, this is important because s8 depends on the mesh size (see Eq. (27)). The next

CR IP T

subsection gives details of the convergence test and the order of accuracy using mesh refinement with the fixed s3 and s8 parameters. The convergence is discussed by comparison with a high order finite volume numerical scheme based on Riemann solver [30, 33].

2.3. Order of accuracy and convergence of the scheme using grid refinement

AN US

Although the order of accuracy of classical D2Q9 [28] is known to be of order

2 in momentum and of order 1 in pressure (density), this section provides a validation of such order through convergence test. Unfortunately, to the authors’ knowledge, there is no available analytical solution for the the multi-dimensional water-hammer case. Therefore, the validation of the order of accuracy is conducted by mesh refinement taking a ”converged” solution as reference (assumed

M

”exact”). The reference solution is considered for the cases where the mesh is refined to Ny = 320 and Nx = 800000. This reference solution is chosen because

ED

the error in this case becomes very small (at order 10−7 ) and the convergence slope remains consistent with the expected order of accuracy. Therefore, the computational mesh is refined where Ny is increased from 20 to 320 and Nx

PT

is increased from 50000 to 800000, respectively (maintaining ∆x = ∆y). The 1 `1 and `2 norms are used to obtain the following errors Em = kpm − pref k1

2 and Em = kpm − pref k2 , where m = {20, 40, 80, 160} is the mesh size, pm is

CE

the pressure solution, and pref is a reference (assumed exact) solution which is for the case where Ny = 320 and Nx = 800000. Figure 12 gives the errors in

AC

logarithmic scale and shows that the slopes of error are 1.2457 of `1 norm and 1.0012 for `2 norm which is consistent with the expected order of accuracy of the used LB scheme. In classical water hammer application, the transient wave is generated by a sudden closure. The pressure rise due to this sudden closure is represented by a vertical line (pressure jump). Such jump can only be (at most) at first order 16

ACCEPTED MANUSCRIPT

of accuracy even if the scheme is at much higher order. Therefore, for classical water hammer test case, the convergence validation through mesh refinement is done i) using data that excludes the pressure jump (therefore testing the order

CR IP T

of the scheme as shown in Table 1 and Figure 12); ii) using data that includes the pressure jump to see how the scheme catches such jump at order 1 (as shown in Table 2). For example, if a high order FV scheme is used, then a slope lim-

iter must be adopted to reduce the scheme order to 1 at the pressure jump (or hydraulic jump).

AN US

Table 1: Validation of the order of accuracy of the scheme using mesh refinement; the pressure

data considered is taken up to 0.95L/cs seconds in time which excludes any pressure jump where rm =

log(Em /E20 ) . log(m/20)

`1

`2

40

1.1272

0.9280

80

1.2350

0.9810

160

1.3748

1.0947

M

m

ED

Table 2: Validation of the order of accuracy of the scheme using mesh refinement; the pressure data considered is taken up to 6L/cs seconds in time which includes the pressure jumps where log(Em /E20 ) . log(m/20)

CE

PT

rm =

m

`1

`2

40

0.8331

0.7228

80

1.0914

0.8368

160

1.2740

1.0034

AC

The 2-d LB scheme is compared with a 2-d 5th order finite volume numerical scheme based on an approximate Riemann solver (FVRS) [30, 33]. The results are shown graphically in Figure 13. In fact, Figures 13 (a)-(d) show a comparison between the FVRS with Ny fixed at 20 and the LB scheme with Ny = 160, Ny = 80, Ny = 40 and Ny = 20, respectively. For clarity, different zoomed

17

ACCEPTED MANUSCRIPT

areas of these figures at different time intervals are shown. Figure 14 gives pressure variation up to L/cs seconds and shows that the LB schemes converges as the mesh is refined. This is mainly because the numerical viscosity of the

CR IP T

LB scheme reduces as the mesh becomes refined (see Eq. (27)). The results of radial wave variations in Figure 14 (a) become identical between both schemes

with the refined mesh. Figure 15 shows the radial wave variation in the time interval [4.2 L/cs , 5.2 L/cs ] seconds. Again, convergence of the LB scheme is

observed in the numerical test results with overall good agreement between the two schemes when Ny = 160 for the LB scheme (Figure 15 (a)). However, Figure

AN US

15 (a) shows that besides the slight difference in amplitude between results from

the two schemes, the result from the LB scheme is also slightly shifted to the right. This is a dispersion effect. In fact, Figure 16 shows clearly that the LB scheme suffers from a dispersion effect especially at low mesh size, where the jump location is shifted from its original location (t = 4 L/cs ). Therefore, if the LB scheme is used to model high modes (radial and azimuthal waves

M

modes), which are dispersive in nature, care should be taken to account for the dispersion effect. In this case, a very fine mesh needs to be used to ensure

ED

accuracy in modeling the behaviour of dispersive high frequency wave. The next section discusses a particular test case in which dispersive waves of the first high

PT

radial mode (M1) are excited and propagate along an infinite pipe system. 2.4. Discussion on the performance of 2nd order LB scheme when modeling dispersive high frequency waves

CE

In this section, a specific numerical test case is designed. This test considers

AC

the pipe system shown in Figure 17 with the following initial conditions:    ρ = ρ0 = 1000 kg/m3 ,       P = p0 = ρ0 g H0 = 1000 × 9.81 × 10 = 9.81 × 104 P a,   (29) Vy = 0, Vx = 0,      D = 0.4 m, L = 100 m,      c = 1000 m/s. s 18

ACCEPTED MANUSCRIPT

The pipe is assumed infinite, and the downstream boundary is located sufficiently far away so that no reflections will interfere with the measured signals. The tests utilize a source located at the upstream boundary (x = 0). An acoustic

CR IP T

wave is generated from the source and only the right-going wave is considered. The generated waveform at the source is shown in Figure 18. This form is cho-

sen because it is smooth and allows the modeler to select the desired frequency

AN US

bandwidth (FBW). Its mathematical form is:   !2  " #  e e  w β β  c  Pf (t)  sin wc (t − = Ps exp −4 log(10) t − ) wc wc βe2    e  where 0 ≤ t ≤ twave = wβc , (30)

where wc = 2π fc is the angular central frequency (in rad/s) with fc being the central frequency (in Hz); Pf is the pressure at the source; Ps = 0.1 p0 is the maximum pressure at the source with p0 being the initial pressure in the pipe; twave is the duration of the wave generated at the source; and βe is a coefficient

M

that controls the FBW. In this work βe = 80 π. The source is considered circular in shape with a given diameter Ds and is located at the pipe centerline. Initially

      P (r)      

PT

inlet) are

ED

at time t = 0 s, the fluid is at rest. The boundary conditions at the source (the   p0 + Pf for 0 < r < Ds , 2 =  p otherwise. 0

CE

(31) Vy = 0,          ∂Vx   = 0. ∂x Pipe wall boundary conditions are similar to the previous 2-d WH test. Because

AC

of the dispersive effect of the LB scheme, the LB mesh is refined so that Ny × Nx = 500 × 62500. This mesh refinement should ensure that the scheme is free from dispersion effects. Parallel computation for the LB scheme is conducted on 30 cores to reduce computational time. Figures 19 (a) and 19 (b) show the pressure variations (x = 25 m, y = 0) obtained from the FVRS and LB

19

ACCEPTED MANUSCRIPT

computational schemes, respectively. The two waveforms observed in Figure 19 correspond to the separated plane wave mode and the first radial mode [30, 31]. Figure 19 shows that, although using a highly refined mesh, the LB scheme still

CR IP T

suffers from large dissipative effects. In comparison with the 5th order FVRS, the LB scheme is inefficient at 2nd order accuracy. In fact, even with the

advantage of LB scheme having optimized parallelization feature, the 5th order FVRS running on a single core requires less computational time (9h4min versus

10h20min for the LB scheme). This demonstrates the importance of using the high order accuracy of the scheme when modeling high frequency (dispersive)

AN US

waves. The effect of high order terms in LB formulations would need further investigation to determine whether better performance might be achieved when modeling HFWs. This work will be discussed in subsequent paper. 2.5. Conclusion

The lattice Boltzmann method with multiple relaxation times is used to

M

model 1-d and 2-d water hammer test cases. For the 1-d case, it is shown that LBM gives the exact solution of pressure variation for inviscid flow. This is

ED

important since, in the literature and common usage, LBM are typically employed to solve only viscous flow and become unstable for inviscid flow test

PT

cases. In fact, the inviscid flow problem can be handled by choosing the tuning √ parameter (α) to be 1, and it is found that this parameter ( α) represents the CFL number. For the 2-d case, the LB scheme is compared with a 5th order FV scheme based on Riemann solver. First the classical water hammer test

CE

is considered with laminar flow. It is found that the relaxation parameter s3 , which relates to the bulk viscosity, governs the dissipation mechanism of radial

AC

waves. On the other hand, the relaxation parameter s8 , which relates to the shear viscosity, governs the dissipation mechanism of axial waves. Overall, the results show that the LB scheme is capable of modeling water hammer applications with good accuracy and performance for low frequency cases. Second, a specific high frequency test case is considered where high frequency smooth waves are injected into an infinite single pipe with a narrow frequency band20

ACCEPTED MANUSCRIPT

width such that the plane wave mode (M0) and the first radial mode (M1) are excited and become separated as they propagate. The results of this numerical investigation show that the LB scheme becomes unstable and its performance

CR IP T

when modeling high frequency waves is significantly impaired. It is presumed that deterioration in stability and accuracy may be due to neglecting high order

terms of the equivalent macroscopic equations considered during the LB numer-

ical formulation. A more detailed analysis of the dissipation mechanism of the LB scheme when modeling high frequency acoustic waves will be discussed in a

AN US

subsequent paper.

Acknowledgments

The authors thank Pierre Lallemand (Beijing Computational Science Research Center, Beijing, China) for helpful discussion during the elaboration of this work.

M

This study is supported by the Hong Kong Research Grant Council (projects

AC

CE

PT

ED

612712 & 612713 & T21-602/15R) and by the Postgraduate Studentship.

21

ACCEPTED MANUSCRIPT

List of figures

AN US

Figure 2: Stencil for the D2Q9 lattice Boltzmann scheme.

CR IP T

Figure 1: Stencil for the D1Q3 lattice Boltzmann scheme.

Figure 3: Reservoir pipe valve system showing the velocity profile.

ED

M

Figure 4: Reservoir pipe valve system.

Figure 5: 1-d LB solutions of the water hammer equation. Dimensionless pressure versus time

PT

at the pipe valve where α = 1 and ν0 = 0.

CE

Figure 6: 1-d LB solutions of the water hammer equation. Dimensionless pressure versus time

AC

at the pipe valve where α = 1/3 and ν0 = 10−6 .

Figure 7: 1-d LB solutions of the water hammer equation. Dimensionless pressure versus time at the pipe valve where α = 1/3 and s = 1.1 (i.e. ν0 = 5000).

22

ACCEPTED MANUSCRIPT

Figure 8: 2-d LB solutions of the water hammer equation. Dimensionless pressure versus time at the pipe valve at the pipe centerline where s3 = 1.75, s8 = 1.4 (i.e. ν0 = 0.49) and

CR IP T

Ny = 100.

Figure 9: 2-d LB solutions of the water hammer equation. Dimensionless pressure versus time

at the pipe valve at the pipe centerline where s3 = 1.1, s8 = 1.99965 (i.e. ν0 = 10−3 ) and

AN US

Ny = 20.

Figure 10: 2-d LB solutions of the water hammer equation. Dimensionless pressure versus time at the pipe valve at the pipe centerline where s3 = 1.45, s8 = 1.99965 (i.e. ν0 = 10−3 )

M

and Ny = 20.

Figure 11: 2-d LB solutions of the water hammer equation. Dimensionless pressure versus

PT

and Ny = 20.

ED

time at the pipe valve at the pipe centerline where s3 = 1.45, s8 = 1.9965 (i.e. ν0 = 10−2 )

AC

CE

Figure 12: Convergence test.

Figure 13: Dimensionless pressure versus time at the valve and at the pipe centerline using Riemann Solver 3rd /5th order (RS) with Ny = 20 and using LBM s3 = 1.45, s8 = 1.99965

(i.e. ν0 = 10−3 ) and different mesh sizes : (a) Ny = 160, (b) Ny = 80, (c) Ny = 40 and (d) Ny = 20.

23

ACCEPTED MANUSCRIPT

CR IP T

Figure 14: Zoomed version of Figure 13 in the time interval [0, L/cs ].

Figure 15: Zoomed version of Figure 13 in the time interval [4L/cs , 6L/cs ].

AN US

Figure 16: Zoomed version of Figure 13 in the time interval [4.2L/cs , 5.2L/cs ].

M

Figure 17: Sketch of unbounded pipe system.

ED

Figure 18: Probing wave from (βe = 40π). (a) time domain ; (b) frequency domain.

PT

Figure 19: Dimensionless pressure variation with time at the pipe centerline (x = 25 m, y = 0) where Ds = 0.2D. (a) obtained using FVRS scheme. (b) obtained using LBM with s3 = 1.45,

AC

CE

s8 = 1.99965.

24

ACCEPTED MANUSCRIPT

Appendix A. LBM scheme D2Q9 scheme is: 1

1

0

λ

−λ

−λ

0

−λ

λ

λ

−λ

−1

−1

2

2

2

−2

−2

1

1

1

2

0

1

−1

−1

0

2

1

−1

−1

1 0

−1

0

0

0

0

1

−1

1

meq 4

=

meq 6

βρ −

=

meq 8

− λy ,

=

M

λ2 ρ

1

−λ

The equilibrium values are given by:    meq = αρ + λ32 ρ (jx2 + jy2 ),   3 jx meq 5 =−λ,   2 2   meq = jx −jy , 7

1

1



  λ    −λ    2    1 .   1    1    0   −1

CR IP T

1

AN US

The moment matrix for the  1 1 1    0 λ 0    0 0 λ    −4 −1 −1   M = 4 −2 −2    0 −2 0    0 0 −2    0 1 −1  0 0 0

3 2 λ2 ρ (jx

j

(A.1)

+ jy2 ), (A.2)

jx jy λ2 ρ .

ED

The equivalent macroscopic equations [14] at order two are: ∂jy ∂ρ ∂jx + + ∂t ∂x ∂y

AC

∂jy ∂t

O(∆t2 ),

(A.3)

PT

∂ jx2 ∂ jx jy α + 4 ∂ρ + + = 6 ∂x ∂x ρ ∂y ρ

+

λ2

=

λ2 ∆t

CE

∂jx ∂t

=

h

−α 6

h

−α 6



1 s3



1 s3

−

1 2



∂ ∂x



∂ ∂y



∂jx ∂x



∂jx ∂x

+

∂jy ∂y

α + 4 ∂ρ ∂ jx jy ∂ jy2 + λ2 + + = 6 ∂y ∂x ρ ∂y ρ = λ2 ∆t

−

1 2

+

∂jy ∂y



+



+

1 3

1 3



1 s8



1 s8

−

1 2

−

1 2



i ∆jx + O(∆t2 ),



i ∆jy + O(∆t2 ),

(A.4)

(A.5)

where ∆f is the Laplacian operator in r space of the function f . The parameter α+4 α is linked to the sound speed cs = λ . The relaxation parameters s3 6 25

ACCEPTED MANUSCRIPT

and s8 are related to the bulk ζ0 and kinematic ν0 viscosities such as ζ0 =     −αλ2 ∆t 1 1 λ2 ∆t 1 1 − and ν = − 0 6 s3 2 3 s8 2 , respectively (e.g. see [28]). To ensure the isotropicity of the dissipation process (ν0 ), s7 and s8 are taken equal.

CR IP T

At second order accuracy, the coefficient β and the relaxation rates s4 , s5 and s6 play no role in the hydrodynamic behavior of the model, however, they are relevant for the stability and the accuracy of the scheme [20].

Remark Using the inverse of moment matrix M (see Eq. (A.1)), the equilibrium

(A.6)

AN US

distribution given by Eq. (A.2) becomes in the f space:    3 9 fieq = ωi ρ + 3V.ci + (V.ci )2 − |V |2 , i = 0, . . . 8, 2 2 where weights ωi are some fixed numbers such that

P8

i=0

ωi = 1 and ρ, V =

(Vx , Vy ) are respectively the density and velocity of the fluid. The truncated equilibrium Eq. (A.6), which approximate the Maxwell distribution, is used because the velocity space is discretized. This approximation comes naturally

M

from Hermite polynomial expansion [41].

ED

References

[1] P. Asinari, T. Ohwada, E. Chiavazzo, A.F. Di Rienzo., Link-wise artificial compressibility method, Journal of Computational Physics, 231, p. 5109-

PT

5143, 2012.

[2] A. Augier, F. Dubois, B. Graille and P. Lallemand, On rotational invariance of Lattice Boltzmann schemes, Computers and Mathematics with Applica-

CE

tions, 67, p 239-255, 2014.

[3] L. Budinski, Application of the LBM with adaptive grid on water hammer

AC

simulation, Journal of Hydroinformatics, 18, 687-701, 2016.

[4] J. M. Buick, C. A. Greated and D. M. Campbell, Lattice BGK simulation of sound waves, Europhysics Letters, 43, p. 235-240, 1998.

[5] M. Bouzidi, M. Firdaous, P. Lallemand, Momentum transfer of a Boltzmann-lattice fluid with boundaries, Physics of Fluids, 13, p. 34523459, 2001. 26

ACCEPTED MANUSCRIPT

[6] M. H. Chaudhry, Applied Hydraulic Transients, 3rd ed., Springer, New York, (2014). [7] S. Chen, G. D. Doolen, Lattice Boltzmann method for fluid flows, Annual

CR IP T

Review of Fluid Mechanics, 30, p. 329-364, 1998. [8] S. S. Chikatamarla, S. Ansumali, I. V. Karlin, Entropic Lattice Boltz-

mann Models for Hydrodynamics in Three Dimensions, Phys. Rev. Lett., 64, p. 010201, 2006.

[9] B. Crouse, D. Freed, G. Balasubramanian, S. Senthooran, P. Lew, L. Mon-

geau, Fundamental Aeroacoustics Capabilities of the Lattice-Boltzmann

AN US

Method, AIAA Paper, 2, p. 006-2571, 2006.

[10] P. J. Dellar, Bulk and shear viscosities in lattice Boltzmann equations, Phys. Rev. E, 64, p. 031203, 2001.

[11] P. J. Dellar, Lattice Kinetic Schemes for Magnetohydrodynamics, Journal of Computational Physics, 179, p. 95-126, 2002.

[12] P. J. Dellar, An interpretation and derivation of the lattice Boltzmann

M

method using Strang splitting, Comput. Math. Applic., 65, p. 129-141, 2013.

ED

[13] H.F. Duan, P.J. Lee, M.S. Ghidaoui, Y.K. Tung, Extended blockage detection in pipelines by using the system frequency response analysis, Journal of Water Resources Planning and Management, ASCE, 138, 55-63, 2011.

PT

[14] F. Dubois, Une introduction au sch´ema de Boltzmann sur r´eseau, ESAIM: Proceedings, 18, p. 181–215, (2007).

CE

[15] F. Dubois, Equivalent partial differential equations of a Boltzmann scheme, Computers and mathematics with applications, 55, p. 1441-1449, 2008.

AC

[16] F. Dubois, Third order equivalent equation of lattice Boltzmann scheme, Discrete and Continuous Dynamical Systems-Series A, 23, number 1/2, p. 221-2482009, a special issue dedicated to Ta-Tsien Li on the occasion of his 70th birthday, doi: 10.3934/dcds.2009.23.221, 2009.

[17] F. Dubois, P. Lallemand, Towards higher order lattice Boltzmann schemes, Journal of Statistical Mechanics: Theory and Experiment, P06006 doi: 10.1088/1742-5468/2009/06/P06006, 2009. 27

ACCEPTED MANUSCRIPT

[18] F. Dubois, P. Lallemand, Quartic Parameters for Acoustic Applications of Lattice Boltzmann Scheme, Computers and mathematics with applications, 61, p. 3404-3416, 2011, doi:10.1016/j.camwa.2011.01.011.

CR IP T

[19] F. Dubois, P. Lallemand, C. Obrecht, M. M. Tekitek, Lattice Boltzmann model approximated with finite difference expressions, Computers & Fluids, doi: 10.1016/j.compfluid.2016.04.013, 2016.

[20] F. Dubois, P. Lallemand and M.M. Tekitek, On a superconvergent lattice Boltzmann boundary scheme, Computers and Mathematics with Applications, 59, p. 2141, (2010).

AN US

[21] F. Dubois, P. Lallemand, M. M. Tekitek, Generalized bounce back boundary condition for the nine velocities two-dimensional lattice Boltzmann scheme, Computers & Fluids, doi.org/10.1016/j.compfluid.2017.07.001, 2017.

[22] M. S. Ghidaoui, On the fundamental equations of water hammer, Urban Water Journal, 1, p.71–83, (2004)

M

[23] M. Geier, M. Sch¨ onherr, A. Pasquali, M. Krafczyk, The cumulant lattice Boltzmann equation in three dimensions: Theory and validation, Comput-

ED

ers and Mathematics with Applications, 70, p. 507, (2015). [24] I. Ginzburg, F. Verhaeghe, D. d’Humi`eres, Two-Relaxation-Time Lattice Boltzmann Scheme: About Parametrization, Velocity, Pressure and Mixed

PT

Boundary Conditions, Communications in Computational Physics, 3, 427478 ,2008.

CE

[25] D. d’Humi`eres, Generalized Lattice-Boltzmann Equations, in Rarefied Gas Dynamics: Theory and Simulations, 159 of AIAA Progress in Astronautics and Astronautics, p. 450-458, 1992.

AC

[26] D. d’Humi`eres, I. Ginzburg, M. Krafczyk, P. Lallemand and L.S. Luo, Multiple-relaxation-time lattice Boltzmann models in three dimensions, Philosophical Transactions of the Royal Society, London, 360, p. 437-451, 2002.

[27] M. Junk, A. Klar, L.S. Luo, Asymptotic analysis of the lattice Boltzmann equation, Journal of Computational Physics, 210, p. 676-704, 2005. 28

ACCEPTED MANUSCRIPT

[28] P. Lallemand, L-S. Luo, Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability, Physical Review E, 61, p. 6546-6562, June 2000.

CR IP T

[29] P. J. Lee, H.F. Duan, M.S. Ghidaoui, B. Karney, Frequency domain analysis

of pipe fluid transient behaviour, Journal of hydraulic research, 51, 609-622, 2013.

[30] M. Louati, In-depth study of plane wave-blockage interaction and analysis of high frequency waves behaviour in water-filled pipe systems (Doctoral dis-

sertation), HKUST, Available at: (http://lbezone.ust.hk/bib/b1618552),

AN US

(2016).

[31] M. Louati, M. S. Ghidaoui, High frequency acoustic waves in pressurized fluid in a conduit. Part 1: Dispersion behaviour, Journal of Hydraulic Research, DOI: 10.1080/00221686.2017.1354931, (2017).

[32] M. Louati, M. S. Ghidaoui, High frequency acoustic waves in pressur-

M

ized fluid in a conduit. Part 2: Range of propagation & preliminary implications in practice, Journal of Hydraulic Research, Submitted. DOI: 10.1080/00221686.2017.1354934, (2017).

ED

[33] M. Louati, M. S. Ghidaoui, The need of high order numerical scheme for modeling dispersive high frequency acoustic waves in water-filled pipe, Journal of Hydraulic Research, under review , (2017).

PT

[34] M. Louati, M. Meniconi, M.S. Ghidaoui, B. Brunone, Experimental study of the eigenfrequency shift mechanism in blocked pipe system, J. Hydraul.

CE

Eng, 10.1061/(ASCE)HY.1943-7900.0001347, 04017044.

[35] S. Mari´e, D. Ricot, P. Sagaut, Comparison between lattice Boltzmann

AC

method and Navier-Stokes high order schemes for computational aeroacoustics, Journal of Computational Physics, 228, p. 1056–1070, 2009.

[36] S. Meniconi, H.F. Duan, P.J. Lee, B. Brunone, M.S. Ghidaoui, M. Ferrante, Experimental investigation of coupled frequency and time-domain transient test-based techniques for partial blockage detection in pipelines, Journal Hydraulic Engineering ASCE, 139, 1033-1040, 2013.

29

ACCEPTED MANUSCRIPT

[37] C. Obrecht, P. Asinari, F. Kuznik, J.J. Roux, High-performance implementations and large-scale validation of the link-wise artificial compressibility method, Journal of Computational Physics, 275, p. 143-153, 2014.

CR IP T

[38] T. Ohwada, P. Asinari, Artificial Compressibility Method Revisited:

Asymptotic Numerical Method for Incompressible Navier-Stokes Equations, Journal of Computational Physics, 229, p. 1698-1723, 2010.

[39] T. Ohwada, P. Asinari, D. Yabusaki, Artificial Compressibility Method and Lattice Boltzmann Method: Similarities and Differences, Computers and

AN US

Mathematics With Applications, 61, p. 3461-3474, 2011.

[40] Y. H. Qian, D. D’Humi`eres, P. Lallemand, Lattice BGK Models for NavierStokes Equation, Europhysics Letters, 17, p. 479-484, 1992. [41] X. W. Shan, X. F. Yuan, H. D. Chen, Kinetic theory representation of hydrodynamics: a way beyond the Navier Stokes equation, J. Fluid Mech. 550, p. 413–441, 2006.

M

[42] J. Wang, D. Wang, P. Lallemand, L-S. Luo, Lattice Boltzmann simulations of thermal convective flows in two dimensions, Computers and Mathematics

ED

with Applications, 65, p. 262-286, 2013. [43] J. Yepez, Quantum Lattice-Gas Model for the Burgers Equation, Journal of Statistical Physics, 107, p. 203-224, 2002.

PT

[44] E. B. Wylie, V. L. Streeter, L. Suo, Fluid transients in systems, Prentice

AC

CE

Hall Englewood Cliffs NJ, (1993).

30

ACCEPTED MANUSCRIPT

AN US

CR IP T

Figure 1: Stencil for the D1Q3 lattice Boltzmann scheme.

AC

CE

PT

ED

M

Figure 2: Stencil for the D2Q9 lattice Boltzmann scheme.

Figure 3: Reservoir pipe valve system showing the velocity profile.

CR IP T

ACCEPTED MANUSCRIPT

CE

PT

ED

M

AN US

Figure 4: Reservoir pipe valve system.

AC

Figure 5: 1-d LB solutions of the water hammer equation. Dimensionless pressure versus time at the pipe valve where   1 and  0  0 .

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 6: 1-d LB solutions of the water hammer equation. Dimensionless pressure versus time at the pipe valve where   1 3 and  0  106 .

Figure 7: 1-d LB solutions of the water hammer equation. Dimensionless pressure versus time at the pipe valve where   1 and s  1.1 (i.e.  0  5000 ).

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 8: 2-d LB solutions of the water hammer equation. Dimensionless pressure versus time at the pipe valve at the pipe centerline where s3  1.75 , s8  1.4 (i.e.  0  0.49 ) and Ny = 100.

Figure 9: 2-d LB solutions of the water hammer equation. Dimensionless pressure versus time at the pipe valve at the pipe centerline where s3  1.1 , s8  1.99965 (i.e.  0  103 ) and Ny = 20.

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 10: 2-d LB solutions of the water hammer equation. Dimensionless pressure versus time at the pipe valve at the pipe centerline where s3  1.45 , s8  1.99965 (i.e.  0  103 ) and Ny = 20.

Figure 11: 2-d LB solutions of the water hammer equation. Dimensionless pressure versus time at the pipe valve at the pipe centerline where s3  1.45 , s8  1.9965 (i.e.  0  102 ) and Ny = 20.

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 12: Convergence test.

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 13: Dimensionless pressure versus time at the valve and at the pipe centerline using Riemann Solver 3rd/5th order (RS) with Ny = 20 and using LBM s3  1.45 , s8  1.99965 (i.e.  0  103 ) and different mesh sizes: (a) Ny = 160, (b) Ny = 80, (c) Ny = 40 and (d) Ny = 20.

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 14: Zoomed version of Figure 13 in the time interval [0, L/cs].

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 15: Zoomed version of Figure 13 in the time interval [4L/cs, 6L/cs].

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

Figure 16: Zoomed version of Figure 13 in the time interval [4.2L/cs, 5.2L/cs].

ED

M

AN US

Figure 17: Sketch of unbounded pipe system.

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

Figure 18: Probing waveform (   40 ). (a): time domain ; (b): frequency domain.

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure 19: Dimensionless pressure versus time at the pipe centerline (x = 25 m away from the source; y = 0) where Ds = 0:2D. (a) obtained using FVRS scheme. (b) obtained using LBM with s3  1.45 , s8  1.99965 .