Parallel Accelerated Fifth-Order WENO Scheme-Based Pipeline Transient Flow Solution Model

Featured Application

Many large-scale water transmission projects, such as China’s South-to-North Water Transfer Project, are hundreds of kilometers long and need real-time regulation applications for management and control. This requires the calculation model to have higher calculation accuracy and faster calculation speed. The existing calculation models are difficult to take these two aspects into account at the same time. In this paper, the accuracy and efficiency of one-dimensional model algorithms are studied and compared, and the best calculation scheme is summarized.

Abstract

The water hammer phenomenon is the main problem in long-distance pipeline networks. The MOC (Method of characteristics) and finite difference methods lead to severe constraints on the mesh and Courant number, while the finite volume method of the second-order Godunov scheme has limited intermittent capture capability. These methods will produce severe numerical dissipation, affecting the computational efficiency at low Courant numbers. Based on the lax-Friedrichs flux splitting method, combined with the upstream and downstream virtual grid boundary conditions, this paper uses the high-precision fifth-order WENO scheme to reconstruct the interface flux and establishes a finite volume numerical model for solving the transient flow in the pipeline. The model adopts the GPU parallel acceleration technology to improve the program’s computational efficiency. The results show that the model maintains the excellent performance of intermittent excitation capture without spurious oscillations even at a low Courant number. Simultaneously, the model has a high degree of flexibility in meshing due to the high insensitivity to the Courant number. The number of grids in the model can be significantly reduced and higher computational efficiency can be obtained compared with MOC and the second-order Godunov scheme. Furthermore, this paper analyzes the acceleration effect in different grids. Accordingly, the acceleration effect of the GPU technique increases significantly with the increase in the number of computational grids. This model can support efficient and accurate fast simulation and prediction of non-constant transient processes in long-distance water pipeline systems.

1. Introduction

The water hammer pressure wave generated when a pipeline suddenly closes a valve can cause pipeline vibration, damage the valve, or cause a pipeline burst in severe cases, resulting in significant engineering accidents. Simulating, predicting, and controlling the water hammer phenomenon are critical research topics in relevant fields [1,2,3]. The commonly used solution schemes, including the MOC and the finite difference method, attempt to calculate water hammers [4,5]. The MOC requires a regular grid, while the spatial step and time step are strictly constrained by the characteristic line, limiting its applicability. Many scholars have improved the classical difference solution schemes. For example, Chaudhry and Hussaini employed the MacCormack scheme to construct an explicit finite difference scheme for solving the pipeline water hammer problem [6]. They found that the second-order finite difference scheme can achieve a computational accuracy higher than the MOC. Wan and others applied the MacCormack scheme of variable grids to transient flow simulation [7]. The mentioned scheme adopts the two-step prediction correction method to realize the flexible division of the spatial grid. Wylie and Streeter proposed an implicit finite-difference scheme for solving the water hammer problem in pipe network systems [8]. Although the main advantage of this scheme is that it remains stable for significant time steps, it requires a time-consuming iteration process. Harten proposed a scheme for solving hyperbolic partial differential equations and proved that this scheme has high resolution and no spurious oscillations at the discontinuities [4]. The finite volume method has been widely utilized in recent years for solving water hammer waves in pipelines [9]. It can provide better conservation performance and less computational time than the TVD (Total Variation Diminishing) difference scheme. Zhao and Ghidaoui constructed the first-order and second-order finite volume methods based on the Godunov scheme using the MUSCL-Hancock method to solve the water hammer problem [10]. They indicated that the second-order accuracy has better conservation at a low Courant number. Waagan proposed the prediction-correction MUSCL-Hancock finite volume method, which guarantees the accuracy of the scheme while ensuring good stability [11].

However, the second-order scheme is still limited in capturing accuracy for this aspect of strongly excited intermittent waves, and better computational results can be obtained using a higher-order scheme. The high-order ENO and WENO schemes significantly improve the capability for intermittent wave capturing [12]. The ENO scheme selects only one template among the multiple templates as the computational template, leading to inefficient computation [13]. In contrast, the WENO scheme utilizes all the templates by assigning corresponding weights to different templates through the convex combination between different templates to improve the intermittent capture performance without causing the waste of templates [14,15,16,17]. For example, Li et al. [18] obtained a new global smoothness index through the Taylor expansion of the local smoothness index. In the framework of the traditional WENO-Z format, an improved third-order finite-difference WENO format was proposed. Li et al. [19] established an improved third-order finite-difference WENO scheme based on convex combinations of different polynomials to calculate the numerical fluxes at the cell boundaries and its truncation errors and parametrization are smaller than those of other third-order WENO schemes. Moreover, the numerical results indicated that the proposed scheme performs better than other third-order WENO schemes.

To judge the performance of the model, in addition to the calculation accuracy of the model, the calculation efficiency of the model is also crucial. If the computational efficiency of the numerical scheme is too low, it will affect the scope of application, such as complex pipe network systems or timely warning of emergency events. Parallel computing is the primary strategy developed in recent years to improve computational efficiency, while the graphics processing unit (GPU) acceleration technology is the direction of rapid development [20,21,22,23]. Liang et al. [24] employed the GPU acceleration technology to establish a high-efficiency and high-precision hydrodynamic model of a rapid movement process of surface water flow for simulating rain and flood in a large-scale watershed resulting in a good acceleration effect. Meng et al. [25] combined the GPU technology with the MOC of a one-dimensional open channel and indicated that the GPU technology could significantly improve the computational efficiency when the number of grids is large. Some scholars applied the GPU technique to solve the system of shallow water equations by taking the WENO scheme and achieved a significant acceleration effect [26,27].

In this paper, the fifth-order WENO finite volume scheme is applied to a one-dimensional pipeline transient flow system, and a pipeline water hammer wave solution model is established based on GPU acceleration technology. This can provide technical support for the intelligent management of complex pipeline systems and rapid prediction of accident emergency dispatch.

2. Control Equations

To calculate water hammers in pipes, the following momentum and continuity equations can be employed [1].

$$g\frac{\partial h}{\partial x}+\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x}+\frac{f}{2D}v\left|v\right|-g\sin \theta =0$$

(1)

$$\frac{\partial h}{\partial t}+v\frac{\partial h}{\partial x}+\frac{c^2}{g}\frac{\partial v}{\partial x}=0$$

(2)

where $h$ is the head of the pressure head, $v$ is the average flow velocity of the pipe section, $f=8g{n}^2{(4/D)}^{1/3}$ is the resistance coefficient along the pipe, $n$ is the roughness, $D$ is the pipe diameter, $x$ is the distance along the pipe axis, and c is the wave speed in the pipe.

Equations (1) and (2) can be expressed in the following linear form:

$$\frac{\partial \mathit{u}}{\partial t}+\mathit{A}\frac{\partial \mathit{u}}{\partial x}=\mathit{s}$$

(3)

where

$$\mathit{u}=\left[\begin{array}{c}h\\ v\end{array}\right]\mathit{A}\mathbf{=}\left[\begin{array}{cc}v& c^2/g\\ g& v\end{array}\right]\mathit{s}=\left[\begin{array}{c}0\\ g\sin \theta -\frac{f\left|v\right|v}{2D}\end{array}\right]$$

Equation (3) can be written in the following approximate form, which can be solved by solving the Riemann problem [28].

$$\frac{\partial \mathit{u}}{\partial t}+\frac{\partial \mathit{F}(\mathit{u})}{\partial x}=\mathit{s}(\mathit{u})$$

(4)

where $\mathit{F}{(\mathit{u})}_x=\overline{\mathit{A}}\mathit{u},\overline{\mathit{A}}=\left[\begin{array}{cc}\overline{v}& c^2/g\\ g& \overline{v}\end{array}\right]$, $\overline{v}$ is the average flow velocity. If the effect of convective terms is ignored, it can be taken as $\overline{v}=0$. This paper takes the arithmetic average proposed by Toro [29] as $\overline{v}=\left(v_i+v_{i+1}\right)/2$. The eigenvalues of the coefficient matrix $\overline{\mathit{A}}$ are ${\lambda }^{\pm }=\overline{v}\pm c$.

Figure 1 shows the schematic diagram of pipeline meshing, and the finite volume method is utilized to attain the discrete solution. ${\mathit{F}}_{i-1/2}$ is the value of $i-1/2$ interface flux and ${\mathit{U}}_i$ is the mean value of $\mathit{u}$ in cell $i$. Taking the control cell $i$ as an example, Equation (4) can be integrated along the x-direction from the control cell interface $i-1/2$ to the interface $i+1/2$ (The range of i is 0 to N, and N is the number of pipe grids).

$$\frac{\partial }{\partial t}{\displaystyle {\int }_{i-1/2}^{i+1/2}\mathit{u}}dx+{\mathit{F}}_{i+1/2}-{\mathit{F}}_{i-1/2}={\displaystyle {\int }_{i-1/2}^{i+1/2}\mathit{s}}dx$$

(5)

Equation (5) is the conservation form of the momentum and continuity equations in the control unit $i$. Let ${\mathit{U}}_i=\frac{1}{dx}{\displaystyle {\int }_{i-1/2}^{i+1/2}\mathit{u}}dx$; if both sides of the Equation (5) are divided by $\Delta x$, the following equation can be obtained (n is a known time step and n + 1 is an unknown time step).

$${\mathit{U}}_i^{n+1}=U_i^n+\frac{\Delta t}{\Delta x}{\displaystyle {\int }_{i-1/2}^{i+1/2}\mathit{s}}dx-\frac{\Delta t}{\Delta x}({\mathit{F}}_{i+1/2}-{\mathit{F}}_{i-1/2})$$

(6)

3. Numerical Calculation Method

This paper takes the WENO-JS format proposed by Jiang and Shu as an example to illustrate how the WENO scheme is applied in a one-dimensional pipeline [14,15,16].

3.1. Numerical Flux Decomposition

The Lax-Friedrichs flux splitting method is employed to obtain the cell boundary flux values as [30]

$${\mathit{F}}_{i+1/2}^{n+1}({\mathit{U}}_i)={\mathit{F}}^{+}({\mathit{U}}_i)+{\mathit{F}}^{-}({\mathit{U}}_i)=\frac{1}{2}({\mathit{F}}^n({\mathit{U}}_i)+\alpha {\mathit{U}}_i)+\frac{1}{2}({\mathit{F}}^n({\mathit{U}}_{i+1})-\alpha {\mathit{U}}_{i+1})$$

(7)

where $\alpha =\max (|\mathit{F}\prime ({\mathit{U}}_i)|)=\max (\left|\overline{v}\pm c\right|)$ is the maximum eigenvalue of the Jacobi matrix $\mathit{A}=\partial \mathit{F}/\partial \mathit{u}$, and it is recommended to take $\overline{v}+c$. The fifth-order accuracy WENO scheme is employed to obtain the flux solution by reconstructing the positive and negative flux vectors ${\mathit{F}}^{\pm }({\mathit{U}}_{\mathit{i}})$.

3.2. Fifth-Order WENO Scheme Flux Reconstruction

According to Equation (7), the numerical flux can be obtained by summing the positive and negative parts, while the following template is chosen for the positive flux ${\mathit{F}}^{+}({\mathit{U}}_i)$: $S_0=\left\{x_{i-2},x_{i-1},x_i\right\},S_1=\left\{x_{i-1},x_i,x_{i+1}\right\},S_2=\left\{x_i,x_{i+1},x_{i+2}\right\}$. Figure 2 shows the location of the templates.

From the above template, the flux reconstruction method of the fifth-order WENO scheme can be obtained as

$${\mathit{f}}_0^n=\frac{1}{3}{\mathit{F}}_{i-2}^n-\frac{7}{6}{\mathit{F}}_{i-1}^n+\frac{11}{6}{\mathit{F}}_i^n$$

(8)

$${\mathit{f}}_1^n=-\frac{1}{6}{\mathit{F}}_{i-2}^n+\frac{5}{6}{\mathit{F}}_{i-1}^n+\frac{1}{3}{\mathit{F}}_i^n$$

(9)

$${\mathit{f}}_2^n=\frac{1}{3}{\mathit{F}}_{i-2}^n+\frac{5}{6}{\mathit{F}}_{i-1}^n-\frac{1}{6}{\mathit{F}}_i^n$$

(10)

The above three templates are employed for weighted reconstruction in the actual calculation, and the weights can be determined using the template smoothness function. This paper employs the template smoothness function $I{S}_k^{\pm }$ proposed by Vukovic and Sopta [15].

$$I{S}_k^{+}={\displaystyle \sum _{\alpha =1}^3{\displaystyle \underset{x_{i-1/2}}{\overset{x_{i+1/2}}{\int }}{\left(\Delta x\right)}^{2\alpha -1}{\left(\frac{d^{\alpha }{\mathit{f}}_k(x)}{d{x}^{\alpha }}\right)}^2}}dx,(k=0,1,2)$$

(11)

Bringing ${\mathit{f}}_k^n$ from Equations (8)–(10) into Equation (11) gives

$$I{S}_0^{+}=\frac{13}{12}{\left({\mathit{F}}_{i-2}^n-2{\mathit{F}}_{i-1}^n+{\mathit{F}}_i^n\right)}^2+\frac{1}{4}{\left({\mathit{F}}_{i-2}^n-4{\mathit{F}}_{i-1}^n+3{\mathit{F}}_i^n\right)}^2$$

(12)

$$I{S}_1^{+}=\frac{13}{12}{\left({\mathit{F}}_{i-1}^n-2{\mathit{F}}_i^n+{\mathit{F}}_{i+1}^n\right)}^2+\frac{1}{4}{\left({\mathit{F}}_{i-1}^n+0{\mathit{F}}_i^n-{\mathit{F}}_{i+1}^n\right)}^2$$

(13)

$$I{S}_2^{+}=\frac{13}{12}{\left({\mathit{F}}_i^n-2{\mathit{F}}_{i+1}^n+{\mathit{F}}_{i+2}^n\right)}^2+\frac{1}{4}{\left(3{\mathit{F}}_i^n-4{\mathit{F}}_{i+1}^n+{\mathit{F}}_{i+2}^n\right)}^2$$

(14)

where the weighting factor ${\omega }_k$ is defined as

$${\omega }_k=\frac{{\beta }_k}{{\displaystyle \sum _{i=1}^3{\beta }_k}}$$

(15)

where ${\beta }_k=\frac{C_k}{{(\epsilon +I{S}_k^{+})}^2},(k=0,1,2)$, and linear combination coefficients are chosen as

$$C_0=0.3,C_1=0.6,C_2=0.1$$

(16)

Based on the template smoothness function, the positive vector flux expression can be derived as

$${\mathit{F}}_{i+1/2}^{+}({\mathit{U}}_i)={\displaystyle \sum _{k=0}^2{\omega }_{k}I{S}_k^{+}({\mathit{F}}_{i-2}^n,{\mathit{F}}_{i-1}^n,{\mathit{F}}_i^n,{\mathit{F}}_{i+1}^n,{\mathit{F}}_{i+2}^n)}$$

(17)

Similarly, the negative flux expression ${\mathit{F}}_{i+1/2}^{-}({\mathit{U}}_i)$ can be obtained from the symmetric distribution of the grid interface $i+1/2$, while the value of the interface flux ${\mathit{F}}_{i+1/2}^{n+1}({\mathit{U}}_i)$ at the moment of $n+1$ can be obtained from Equation (7).

3.3. Time Layer Discrete Method

The above steps can only find the value of the interface flux at $n+1$ times, while the value of the source term $\mathit{s}(\mathit{u})$ is required to calculate the ${\mathit{U}}_{{}_i}^{n+1}$ at $n+1$ times. Now, the second-order Runge-Kutta method is employed to solve the initial value problem.

First, the flow term is updated as

$${\overline{\mathit{U}}}_i^{n+1}={\mathit{U}}_i^n+\frac{\Delta t}{\Delta x}(\mathit{F}({\mathit{U}}_{iL})-\mathit{F}({\mathit{U}}_{iR}))$$

(18)

In the second step, the source item time step is updated as

$${\mathit{U}}_i^{n+1}={\overline{\mathit{U}}}_i^{n+1}+\Delta t\mathit{s}({\overline{\mathit{U}}}_i^{n+1}+\frac{\Delta t}{2}\mathit{s}({\overline{\mathit{U}}}_i^{n+1}))$$

(19)

3.4. Boundary Condition Processing

In this paper, the treatment methods of upstream and downstream boundary conditions and pipeline bifurcation boundary are given.

3.4.1. Upstream Boundary Conditions

In actual engineering, the upstream boundary of the pipeline is primarily constant pressure boundary conditions, such as pressure front pool and constant reservoir level. Combining the boundary conditions with MOC gives

$$\left\{\begin{array}{l}{v}_{(1,1)}^n=\frac{g}{c_{}}\left(h_{(1,1)}^n-h_{(1,0)}\right)+v_{(1,0)}(1+\frac{g\Delta t\sin \theta }{c_{}}-\frac{f\Delta t}{2D_1}\left|v_{(1,0)}\right|)\\ h_{(1,0)}=h_0\end{array}\right.$$

(20)

where $h_0$ is the upstream constant water pressure and $c_{}$ is the pipeline pressure wave speed.

3.4.2. Downstream Boundary Conditions

The downstream boundary can be the flow rate boundary, flow boundary, or valve boundary. The downstream flow rate is given as $v_{(i,N(i))}$ or flow rate as $\left(v_{(i,N(i))}=Q_{(i,N(i))}/A_i\right)$ for the flow rate boundary or flow boundary (Here i represents the pipe number of the pipe network, the same as below). Now, the head at the end of the pipe can be calculated as

$$h_{(\mathrm{i},N(i))}=h_{(\mathrm{i},N(i)-1)}^{n-1}-v_{(i,N(i))}\frac{c_i}{g}+v_{(\mathrm{i},N(i)-1)}^{n-1}\left(\frac{c_i}{g}+\Delta t\sin \theta -\frac{f\Delta x\left|v_{(\mathrm{i},N(i)-1)}^{n-1}\right|}{2g{D}_i}\right)$$

(21)

Based on the MOC at the boundary and the orifice outflow law, the relationship between head and flow velocity under a valve boundary condition can be determined as

$$v_{(i,N(i))}=\frac{v_m}{\sqrt{h_0^{\prime }}}{\tau }_j\sqrt{h_{(\mathrm{i},N(i))}}$$

(22)

where $v_m$ stands for the valve fully open valve at the average flow rate, $h_0^{\prime }$ stands for the pipeline constant flow state at the valve head, and ${\tau }_j$ is the valve opening degree.

To obtain the upstream and downstream boundary flux values ${\mathit{F}}_{-1/2}$ and ${\mathit{F}}_{N+1/2}$, the virtual boundary values outside the cell inside the pipe should be determined. The dummy grid technique proposed by Mohamed et al. [31] is employed in this paper to improve the computational efficiency while ensuring the accuracy of the results. As shown in Figure 1, for the upstream constant pressure boundary, we have ${\mathit{U}}_{-2}={\mathit{U}}_{-1}={\mathit{F}}_{-1/2}$, while for the downstream flow velocity boundary or valve boundary, we have ${\mathit{U}}_{N+2}={\mathit{U}}_{N+1}={\mathit{F}}_{N+1/2}$.

3.4.3. Bifurcated Pipe Interface Handling

When there is a bifurcation in the water delivery system, as shown in Figure 3, it can be handled as the following.

According to the positive characteristic line of pipe A and the negative characteristic line of pipes B and C, we have

$$\left\{\begin{array}{l}{h}_{(A,N(A))}+\frac{c_A}{g}{v}_{(A,N(A))}=h_{(A,N(A)-1)}^{n-1}+v_{(A,N(A)-1)}^{n-1}\left(\frac{c_A}{g}+\Delta t\sin {\theta }_A-\frac{f\Delta x\left|v_{(A,N(A)-1)}^{n-1}\right|}{2g{D}_A}\right)\\ h_{(B,0)}-\frac{c_B}{g}{v}_{(B,0)}=h_{(B,1)}^{n-1}-v_{(B,1)}^{n-1}\left(\frac{c_B}{g}+\Delta t\sin {\theta }_B-\frac{f\Delta x\left|v_{(B,1)}^{n-1}\right|}{2g{D}_B}\right)\\ h_{(C,0)}-\frac{c_C}{g}{v}_{(C,0)}=h_{(C,1)}^{n-1}-v_{(C,1)}^{n-1}\left(\frac{c_C}{g}+\Delta t\sin {\theta }_C-\frac{f\Delta x\left|v_{(C,1)}^{n-1}\right|}{2g{D}_C}\right)\end{array}\right.$$

(23)

The following relation can be obtained from the conservation of flow at the bifurcation and the equal head relationship.

$$A_A{v}_{(A,N(A))}=A_B{v}_{(B,0)}+A_C{v}_{(C,0)}$$

(24)

$$h_{(\mathrm{A},N(A))}=h_{(\mathrm{B},0)}=h_{(C,0)}$$

(25)

Equations (23)–(25) should be solved to derive the flow rate and head value at the bifurcation connection of the corresponding moment. If there are more than three sections of pipes connected at the bifurcation, referring to this paper, the number of pipes should be expanded for calculation.

4. GPU Acceleration Implementation Method

Since the pipeline length can be hundreds of kilometers for large water transmission projects, traditional CPU computing architecture is time-consuming and inefficient. With the introduction of GPU acceleration technology, the model calculation efficiency can be significantly improved with minimal loss of model accuracy, ensuring both calculation accuracy and efficiency.

4.1. CPU and GPU Related Parameters

Due to the rapid development of science and technology, the CPU and GPU are changing year by year, while different models of GPU graphics cards have different acceleration effects. The CPU and GPU used in the model are new products in recent years. The CPU is the latest Core series of 10th-generation machines, and the GPU contains 1536 CUDA cores.

Table 1. Key parameters of GPU and CPU.

GPU Type	Computational Framework	Number of Transistors	Number of Stream Processors	Memory Capacity	Single-Precision Floating-Point	Video Memory Bandwidth
NVIDIA GeForce GTX 1660Ti	Pascal	6 billion 600 million	1536	6000 MB	4.85 TERAFLOPS	288(Gb/s)
CPU Type	Computational Framework	Basic Frequency	Number of Threads	L3 Cache	Maximum Memory Support
Intel(R) Core (TM) i7-10700	Comet Lake-S	2.90 GHz	8 threads	16 MB	8 GB

shows the specific key parameters of the CPU and GPU.

4.2. GPU Accelerated Parallel Computing Process

The model adopts C++ and CUDA language programming to implement the GPU parallel acceleration process. The initial data is stored and initialized on the CPU side, and the cudaMemcpy function is utilized to transfer the computational data to the GPU side as fast as possible, and the GPU threads are called through the kernel function for parallel processing. Figure 4 depicts the specific computation flow.

5. Model Validation and Case Analysis

To evaluate the model’s computational accuracy and efficiency, this paper selects typical working conditions for simulation and validation. It selects the MOC and the first-order and second-order Godunov schemes [10] for comparison and analysis. The Courant number sensitivity of the model is qualitatively verified through the dissipation analysis of different schemes under the same Courant number. By increasing the number of grids, the acceleration effect of the GPU technology in dividing different grids is studied.

5.1. Parameter Sensitivity Analysis

By calculating and evaluating the numerical dissipation and surge intermittent capture performance of the pressure wave curve at the pipe’s valve at different Courant numbers, the sensitivity of different simulation methods to the Courant number is verified.

Assuming that the two calculation conditions are recorded as Case 1 and Case 2, the pipeline length is 1960 m and 39,200 m, the calculation time is 40 s and 800 s, the grid number is 20 and 40, and the Courant number is 0.1 and 0.01, respectively. Other parameters are the same. The pipe diameter is 1.0 m, the initial flow velocity is 0.5 m/s, the upstream constant water level is 10 m, the wave speed is 980 m/s, and the gravitational acceleration is 9.806 m/s². The viscous term is ignored, so the roughness of the pipe is 0, and the end valve is suddenly closed.

Figure 5 shows the pressure fluctuation characteristics at the end valve and the 540 m upstream of the pipeline when the model Courant number is 0.1 and the model grid is 20 (case 1). According to the Joukowsky theorem, the change value of water hammer pressure:

$$\Delta h=\frac{c\Delta v}{g}=\frac{980\times 0.5}{9.806}=49.97\mathrm{m}$$

(26)

The maximum pressure is 59.97 m, the minimum pressure is −39.97 m, and the phase length (the time for traveling one wavelength) of the pressure wave is $T=2l/c=2\times 1960÷980=4\mathrm{s}.$

As shown in Figure 5, the calculation results of the MOC and the first-order Godunov scheme indicate significant numerical dissipation with increasing time, while the peak pressure waves captured by both schemes are decreasing. The extreme value of the pressure wave at the end valve can be captured accurately during the calculation period using the second-order Godunov scheme. The calculation results at 540 m upstream exhibit a significant error, and the scheme’s numerical dissipation is still more severe. Simultaneously, the second-order Godunov scheme cannot capture the abrupt interruptions. In contrast, the fifth-order WENO scheme accurately captures the actual pressure wave peak and provides performance superior to the other three schemes in capturing the abrupt interruptions of the excitation wave. Figure 5a shows that the second-order Godunov scheme captures a complete excitation interruption with a response time of about 4 s, while the fifth-order WENO scheme captures a complete interruption within 1 s. Besides, the non-physical oscillations generated by the fifth-order WENO scheme are so minimal that their impact on the actual results can be neglected.

Figure 6 shows the pressure fluctuation characteristics when the Courant number is taken as 0.01 and the model grid is 40 (Case 2). Similarly, the pressure change value is $\Delta h=c\Delta v/g=49.97\mathrm{m},$ and the phase length of the pressure wave is $T=2l/c$ $=2\times 39200÷980$ $=80\mathrm{s}.$ Compared with Figure 5, the numerical dissipation of the MOC and the first-order Godunov scheme decreases slightly by increasing the number of computational grids. However, the scheme cannot capture the pressure wave extremes after 100 s, and the scheme calculation results are almost distorted. The second-order Godunov scheme can better capture the intermittent excitation waves and slightly reduce the dissipation of the scheme values. The fifth-order WENO scheme has no significant nonphysical oscillations and numerical dissipation while maintaining high performance in capturing intermittent excitation waves. In summary, the fifth-order WENO scheme has lower numerical dissipation and better intermittent capture capability than the traditional MOC and Godunov scheme for calculating water hammer waves in pipes with low Courant number conditions.

5.2. Model Application and GPU Acceleration Performance Evaluation

The idealized long-distance water pipeline network shown in Figure 7 is constructed to evaluate the simulation’s computational accuracy and the GPU program’s acceleration effect. It is assumed that all the pipelines are long-distance transmission pipes of several tens of kilometers, and points A and B are the branch points of the trig point. The relative elevation of the upstream reservoir constant water level Z1 and Z2 and the downstream water level Za and Zb are 85 m and 5 m, respectively, while the regulating valve at the end of the pipeline is closed by linear variation with a closing time of the 60 s. The MOC and second-order finite volume Godunov scheme are employed to validate the model. The calculation time step is taken as 0.01 s, and the gravitational acceleration is chosen as 9. 806 m/s². The pipeline model parameters are shown in

Table 2. The pipeline model parameters.

Pipe Number	Pipe Length/km	Pipe Diameter/m	Roughness	Wave Speed/(m/s)	Initial Flow/(m3/s)
Pipe 1	39.75	0.981	0.013	996	0.47
Pipe 2	39.75	0.981	0.013	996	0.47
Pipe 3	99.40	1.389	0.012	1000	0.94
Pipe 4	59.64	0.981	0.014	994	0.47
Pipe 5	59.64	0.981	0.014	994	0.47

.

5.2.1. Comparative Analysis of the Results of Different Simulation Methods

Figure 8a shows the process of pressure wave curve change at valve C when the end valves C and D are simultaneously closed for the working condition shown in Figure 7. The working condition first runs steadily for 800 s, and then the valve closes linearly within the 60 s. This paper adopts the fifth-order WENO finite-volume method to take 0.01 and 0.005, respectively. As shown in Figure 8a, the calculation results of different schemes are almost the same, and the error of the extreme value point is within 1‰ when the WENO scheme and the MOC take 0.005. This means that the fifth-order WENO scheme for solving the pipeline water hammer wave can achieve a high accuracy under a very low Courant number. Figure 8b compares the solution results after taking GPU acceleration with the CPU program’s results. The GPU is programmed with double-precision floating-point numbers to ensure calculation accuracy. The maximum error between the calculation result and the CPU program calculation result is 0.5‰, significantly below the allowable error in actual engineering.

Figure 9 shows the flow and water hammer wave variations at different locations when the end valves C and D are simultaneously closed. As shown in Figure 9, due to the symmetry of points C and D location and the simultaneous closure of the valve, the pressure wave curves at valves C and D overlap entirely. Simultaneously, after the valve is closed linearly, the water hammer pressure wave reaches its maximum value (131.9 m). The flow rate at A gradually decreases to 0 with the closure of the valve, while the flow rate at B and C still decrease cyclical in the fluctuation changes due to the role of pipe friction, which gradually decreases the amplitude of fluctuations. Figure 10 shows the flow and pressure head variations at different locations when the valve is closed only at C (linear closure within the 60 s). As shown in Figure 10a, after closing valve C, the pressure head at C suddenly rises, and its peak reaches 95 m, then oscillates down and stabilizes at 64.3 m. After about 40 min, due to the horizontal arrangement of the pipe, the flow rate in pipe 4 drops to zero, while the pressure head at point C is equal to that at point B. The maximum pressure head generated by the whole process system is still at the end of valve C. It can be seen that the downstream valve is still the most likely place to generate damage, even for a complex bifurcated piping system.

5.2.2. Analysis of Model Computation Efficiency and GPU Acceleration Performance

The water hammer change process at the end valve was calculated in time steps of 0.02 s, 0.015 s, 0.01 s, and 0.005 s. The model’s computational efficiency can be derived by analyzing the computational time consumed by the different schemes, with a computational time of 4800 s. According to the CFL constraint, the MOC and the first-order Godunov scheme should meet Cr = 1, while the second-order Godunov scheme takes Cr = 0.5 to obtain the best calculation effect [10], and the proposed fifth-order WENO finite volume scheme can meet Cr > 0.01. The total computation times for different schemes are presented in

Table 3. The calculation time of different simulation methods.

Time Step (t/s)	MOC (t/min)	First-Order Godunov (t/min)	Second-Order Godunov (t/min)	Fifth-Order WENO (t/min)
0.020	3.99	7.76	9.55	3.05
0.015	6.72	13.47	16.75	6.34
0.010	16.31	28.73	36.90	14.36
0.005	59.00	118.75	146.54	22.47

and Figure 11. According to Table 3, the second-order Godunov scheme takes the highest time, while the MOC and the first-order Godunov scheme take the middle time, and the WENO scheme takes the least time for different time steps, while the Courant number meets the scheme requirements. The fifth-order WENO scheme takes 35% to 90% of the MOC and 15% to 39% of the time of the second-order Godunov scheme when the calculation accuracy is guaranteed and the appropriate Courant number is chosen. It can be concluded that the fifth-order WENO finite volume scheme proposed for the water hammer solution model can significantly improve computational efficiency.

The computational time consumption of CPU programs and GPU programs at different grids was analyzed to evaluate the acceleration effect of the GPU-WENO scheme. The CPU model is Intel(R) Core (TM) i7-10700, and the GPU graphics card is NVIDIA GeForce GTX 1660Ti, and both of them are the latest electronic products in recent years. Still using the ideal operating conditions shown in Figure 7, the pressure wave variation during sudden closure of the end valve was calculated for a total time of 48 s, as presented in

Table 4. Effect analysis table of GPU acceleration.

Number of Grids	Grid Accuracy/(m)	CPU-WENO (t/s)	GPU-WENO (t/s)				Acceleration Ratio/(Times)
			Number of Threads				Number of Threads
			8	16	32	64	8	16	32	64
600	400	4.44	5.22	5.21	5.19	5.22	0.85	0.85	0.86	0.85
6000	40	48.16	11.55	6.33	5.20	5.23	4.17	7.61	9.26	9.21
6 × 104	4	479.24	91.28	57.88	47.27	46.96	5.25	8.28	10.14	10.21
6 × 105	0.4	4731.08	793.81	502.24	358.96	354.39	5.96	9.42	13.18	13.35

. As shown in Table 4, when the grid number is uniformly divided into 600 and the grid length is 400 m, the single-threaded computation efficiency of the 8-core CPU is higher than that of the GPU. The reason is that the number of CPU program cycles is small, and the CPU core performance can be fully utilized for the small number of grids, while the GPU program is so time-consuming compared to the CPU program because the internal data is copied from the host side to the device side and then copied from the device side to the host side after completing the calculation. The computational efficiency of the GPU program is significantly higher than that of the CPU program when the number of grids reaches thousands and above. At this point, the CPU program’s large number of round-robin calculations is a critical factor in the program’s time consumption, while the efficiency of the GPU is superior to the CPU program by calling multiple threads on the kernel block to perform parallel calculations. Figure 12 shows that when the number of threads reaches 32 or more, the GPU acceleration effect does not change with the increase of threads, indicating that for engineering problems within hundreds of thousands of grids, 32 threads can achieve the best acceleration effect. Besides, the acceleration effects of different graphics cards are significantly different. This model employs the NVIDIA GeForce GTX 1660Ti graphics card, while the acceleration effect can be more than 13 times compared with a single CPU program by using GPU parallel computing technology with a grid number of 500,000 and a grid accuracy of 0.4 m. With the GPU-WENO scheme, which has a high application prospect in solving practical engineering problems, the computation speed can easily reach more than tens of times that of the commonly used MOC or differential scheme.

6. Discussion

This study indicates that the high-precision scheme can well simulate the propagation process of the pressure wave in the one-dimensional pipeline system. Compared to the MOC and the second-order Godunov scheme, the fifth-order WENO scheme yields lower numerical dissipation and no spurious oscillations. GPU acceleration technology is introduced, and the acceleration effect is reached more than ten times. These results can be applied to the regulation and control of long-distance water conveyance systems. This study does not consider the effect of adding gas into the pipeline system. This issue will be further studied in the future.

7. Conclusions

Based on the one-dimensional pipeline non-constant flow equation set, the virtual boundary technique is adopted at the upper and lower ends of the model. The bifurcated pipeline boundary solution method is combined with the Lax-Friedrichs splitting method for flux decomposition, and the time term is discretized by the second-order Runge-Kutta method. Moreover, a fifth-order WENO finite volume scheme water hammer solution model is established for the long-distance water transmission system. The key results are as follows.

(1) The model is verified to ensure low numerical dissipation and lack of spurious oscillations at tiny Courant numbers under ideal conditions. The model is suitable for large and complex water conveyance systems. Compared to the commonly used MOC, the computation load of the model is 35% to 90% less than the MOC and 15% to 39% less than the second-order Godunov scheme while maintaining the accuracy of the results.

(2) Using GPU acceleration technology to accelerate the model, the model can achieve a significant acceleration effect while ensuring accuracy when the number of grids is greater than thousands. As the number of meshes increases, the GPU acceleration effect becomes more and more significant. The study found that GPU with 32 threads can deal with engineering problems within hundreds of thousands of grids with the best acceleration effect.