A lattice Boltzmann method (LBM) is utilized to solve single-phase transient flow in pipes. In order to eliminate grid limitation related to the method of characteristics, governing equations are modified using appropriate coordinate transformation. The introduced modification removes connection between Courant number and spatial disposition of the computational nodes, forming a more flexible and robust mathematical base for numerical simulations. The computational grid is configured independently of the wave speed, significantly decreasing the demand for computational resources and maintaining the required accuracy of the method. Thereafter, the appropriate equilibrium distribution function for the D1Q3 lattice has been defined. In order to give a comprehensive base for modeling transient flow in complex pipeline systems, detailed elaboration of the corresponding boundary conditions has been given. Two benchmark problems with the corresponding error analysis are used to validate the proposed procedure.

## INTRODUCTION

Transition between two steady states in pipeline systems is a common phenomenon in engineering practice. Caused by sudden change in flow regime (pump failure, instantaneous valve closure), transient flow is characterized as highly unsteady flow with intensive pressure and velocity fluctuations, classically known as water hammer. In order to describe this phenomenon, a set of two hyperbolic partial differential equations, mass and momentum equation is utilized (Fox 1977; Wylie & Streeter 1978; Chaudhry 2014). Due to their hyperbolic nature, these equations are transformed to finite differential equations by using the method of characteristics (MC) (Wylie & Streeter 1978). Hartree (1958) was the first to introduce spatial interpolation in the case where the Courant number (Cr) was less than 1.0, while Wiggert & Sundquist (1977) and Vardy (1976) solved the pipeline transients for characteristics projecting over the fundamental grid size. Their analysis shows the effects of interpolation, spacing, and grid size on numerical attenuation and dispersion. Instead of widely used spatial interpolation, Goldberg & Wylie (1983) used interpolations in time, while Sibetheros *et al.* (1991) investigated the MC with spline polynomials for interpolations. Besides the MC, finite-difference (FD) and finite-volume (FV) methods have also been used in order to solve transient equations. Chudhury & Hussaini (1985) solved water hammer equations by MacCormack, Lambda, and Gabutti explicit FD schemes, while Verwey & Yu (1993) applied a space-compact high-order implicit scheme. Furthermore, Zhao & Ghidaoui (2004) formulated, applied, and analyzed first- and second-order explicit FV Godunov-type schemes for water hammer problems, while Leon *et al.* (2007, 2008) and Leon & Oberg (2013) applied FV method on two-phase water hammer flows. Application of the MC is proven to be simple to code, is accurate and efficient. However, major drawbacks related to the MC arise in practical applications. Since the spatial disposition of computational points is limited to relation , where *a* is the propagation (wave) speed and is time step, frequently a single value for is adopted for the whole system, as being representative. This distance is often prescribed by the system configuration; distance between two inner boundary conditions (for example, two air chambers). Consequently, a long pipe section can impose a large number of computational points, which can greatly affect the efficiency of the method. On the other hand, if spatial or time interpolation is used, additional procedures are inevitably introduced, which again influence the stability and accuracy of the method (Vardy 1976; Wiggert & Sundquist 1977; Goldberg & Wylie 1983).

In order to somehow overcome deficiencies related to the previously described methods, a lattice Boltzmann method (LBM) (Wolf-Gladrow 2005) is offered as an alternative approach. The first application of the LBM on transient flow was considered by Cheng *et al.* (1998), while only a case of rather simple practical implementation was done by Wu *et al.* (2008). In this paper, further development of the LBM is considered. In order to enhance the efficiency of the method, adaptive grid procedure is introduced. With this modification, restriction between *Cr* and spatial step is removed, and arbitrary grid configuration approach is established. Hence, the number of computation points now can be significantly reduced, which is a desirable feature for long pipe sections. The corresponding grid size accuracy analysis is also conducted. Furthermore, mathematical formulation of most used boundary conditions is presented in detail. Implementation of the basic pipeline elements to the lattice Boltzmann model, such as pump, air chamber, valve, expansion/contraction of pipes, branching of pipes, is elaborated from the distribution function point of view, and some simplifications are introduced. This will much enhance the applicability of the method, with the opportunity to efficiently simulate complex pipeline systems. In contrast to the few lattice Bolzmann applications on transient flows available in the literature, where only the basics of the LBM are presented, the adaptive grid approach followed by detailed and comprehensive elaboration of the boundary condition and practical implementation offers scientists and engineers a robust and efficient tool for managing the transient flow in complex pipeline systems. The presented method was tested and verified on two different examples, while the MC was used for comparison.

## TRANSIENT FLOW EQUATIONS

*t*is time,

*x*is Cartesian coordinate, is hydraulic head,

*V*is longitudinal mean velocity,

*g*is gravitational acceleration, and is friction factor. Wave speed

*a*is defined as:where

*K*is bulk modulus,

*E*is Young's modulus, is volumetric mass density and

*e*is wall thickness. Depending on the pipe anchors, parameter is calculated as:

– pipe anchored upstream,

– pipe anchored throughout (no axial dilatation),

– expansion joints (axial dilatation allowed),

## LATTICE BOLTZMANN MODEL

*et al.*1954) is utilized. The evolution equation is defined as:where is the particle distribution function along the link, is the local equilibrium distribution function, is the position vector in the 1D domain,

*t*is time, is the time step, is the force term, and is relaxation time. It should be noted that Equation (6) is defined using the new -domain, therefore the required symmetry for discrete particle velocities is ensured. For the three-velocity lattice, particle velocities along the direction are defined as and . Accordingly, is lattice size in direction.

*et al.*(2010), the following formulation is proposed here:Term

*G*in Equation (7) actually denotes connection between the physical and computational (lattice) domains, and it is calculated prior to the LBM computations. Furthermore, discrete formulation takes the form , where is physical distance between computational nodes.

*V*:

### Deduction of transient flow equations

### Boundary conditions

Water supply systems are often characterized by a variety of elements (fittings, pumps, and storage tanks, etc.), hence adequate implementation of the boundary conditions represents one of the major tasks when transient pipe flow modeling is considered. Contrary to open channel flow cases, where formally only two types of boundary conditions are used, in transient pipe flow every type of element requires appropriate and specific LB formulation. Therefore, the detailed presentation of boundary conditions in the case of the transient LBM is presented in the following section.

### Contraction/expansion of a pipe

*Q*between two boundary sections of the pipes and are used (Figure 1):

*j*denotes the corresponding pipe, while index marks the computational node along pipe

*j*. It is obvious from Figure 1 that and are the unknown distribution functions related to pipes and , respectively. Introducing Equation (10) into Equations (21) and (22) and then expressing them in terms of unknown distribution functions and , respectively, giveswhere

*A*is the cross section area of the corresponding pipe. Finally, substituting Equation (24) into Equation (23) leads to final expressions for calculating the unknown distribution functions in the formsIt should be noted that Equations (25) and (26) are equally applicable for both cases, i.e., contraction and expansion, which makes them rather universal when change of pipe diameter is considered.

### Valve

### Valve located at the end of the pipe

### Branching pipes

### Reservoir located at the end of the pipe

*T*, providing constant hydraulic head, and cross section of the pipe is introduced (Figure 5):Applying zeroth moment for hydraulic head and first moment for discharge (Equation (10)), the following equation is derived:It is evident that Equation (41) is quadratic; hence the solution for the unknown distribution function is defined as:whereFurthermore, for a reservoir located at the right end of the pipe, the same procedure is utilized. Using the energy equation in the formthe solution of the corresponding quadratic equation is defined as:where

### Surge tank

### Air chamber

*m*is the polytropic index. Introducing the finite difference form of Equations (56) and (57) along with Equation (59) in to Equation (58), relation for the pressure head inside the chamber is obtained:where

*n*denotes time level. Further, combining Equation (50) – representing relationship between flow in the chamber and hydraulic head in the pipe – and Equation (60) with Equation (55), a quadratic equation similar to Equation (51) is derived:

Finally, solution is obtained using Equation (52) with the following replacements , and .

### Pump

*l*is the moment of inertia of the pump and entrained liquid, and is the change in angular speed with time. It should be noted that in the case of pump failure the torque of the pump becomes zero instantaneously. By introduction of relations and , Equation (66) can be discretized using the finite difference method as:where is derived applying the four quadrant relation in the formThe set of Equations (62), (63), (64), (67), and (68) close the system of required equations, which is then solved iteratively. The procedure is described here:

In the first step, functions and are determined using assumed values for

*v*and .In the second step, dimensionless speed of rotation and torque is calculated using Equations (67) and (68), respectively.

- Introducing the zeroth statistical moment in the form into Equation (63), and using assumed value in the first iteration for the second unknown distribution function , in the third step the following equation is solved
- In the fourth step, the second unknown distribution function is calculated using the combination of Equation (62) and the first statistical moment
Finally, dimensionless flow is calculated using the distribution function from the previous step.

The entire procedure is repeated from stage 2 using tolerance .

## THE NUMERICAL SIMULATION AND THE VALIDATION OF THE MODEL

In order to test the proposed form of the LBM on more practical transient problems, we focus our test analysis mainly on the pump failure cases. Being the most common reason for inducing transient flows in pipeline systems, pump failure is, at the same time, complex enough and representative enough for testing the novel numerical method.

### Single pipe system

^{3}/h, actual developed head m, efficiency , and speed of rotation rpm is set up (Amarex 2014). From the open reservoir, having water level m, water is pumped along the horizontal pipe vertically positioned at level m. Furthermore, two valves are implemented in the system. However, primarily the role of valves in this example is not flow regulation using an adapted closing dynamics, but rather boundary condition implementation on one hand (free outflow at the end of the pipe), and enforcement of the MC to define the computational grid along the pipe according to the smallest distance limitation (valve located m from the pump) on the other hand. As stated in the Introduction, MC utilizes an approach where the governing equations are solved along the characteristics (trajectories), which can connect two computational points (), or it can be projected outside () or inside () the computational cell. In order to avoid additional computational work and reduce numerical error (caused by interpolation) as much as possible, especially when nonlinear definition for characteristic is used, most commercial softwares utilize a standard model, i.e., when . However, this approach also introduces a certain type of limitation if the computational grid is defined according to the smallest distance in the system. This can significantly increase dimension of the computational grid, which in turn influences overall computational efficiency. For the purpose of testing and verification of the transient LB method, a commercial software AFT Impulse 4.0 (Applied Flow Technology 2011) is used in our investigation.

### A small pipeline system

_{S}). Comparison of the obtained results in the form of hydraulic head variations and hydraulic grade-lines related to the extreme values () of hydraulic head is presented in Figures 15 and 16. Similar to the single pipe system, very good agreement between the two models is achieved. For the time variations of the hydraulic head (Figure 15), small deviations are present only at the local minimum, starting at s. Relative to the maximal drop of hydraulic head, this digression falls in the range of . However, for the minimal and maximal extrema, i.e., for extreme pressure drops instantaneously after pump fails ( s in Figure 15), excellent agreement is obtained. Furthermore, this trend is relatively maintained along the whole system, which is shown in Figure 16. Only at a few computational points, located near the end of the system, is some divergence noted. Also, a 33.22% decrease in simulation time when the LB method is used is achieved. However, this model did not include additional procedures and coding for high performance computing (CUDA, OpenCL), which can further accelerate overall procedure metrics.

^{3}. The chamber is connected to the system via a short pipe having loss coefficient , according to Equation (52). The same procedure of transient regime induction and analysis is applied here as in the previous example. Comparison of the results obtained by the LBM and MC in terms of hydraulic head variations and envelopes is given in Figures 17 and 18, respectively. Once again, very good agreement between the compared results is achieved in both cases. Minimal deviations noted in Figure 18, ranging from relative to the maximal head drop, demonstrate that the LBM model is more than capable for solving a variety of complex practical problems.

## CONCLUSIONS

To overcome grid limitation typical to the MC, requiring Cr equal to one for stability in solving transient equations of pressurized flow, the same problem is solved by the LBM. To provide a higher level of flexibility in grid use, an adaptive grid approach is exploited. Sectional autonomy in the application of specific grid density to each pipe section significantly increases the efficiency of the overall procedure, maintaining the basic stability and accuracy. Introducing the LB method as an alternative to the MC in water hammer calculations provides the opportunity for parallel processing, which further increases the efficiency of calculations and cuts down calculation time. This is of high importance in large and complex pipeline systems. A detailed overview and elaboration of internal and external boundary conditions, i.e., pipeline elements and components is given. A real pipeline system with pump failure is modeled for the verification of the proposed LBM. Very good agreement with the MC is achieved.