Lattice Boltzmann simulation of the open channel flow connecting two cascaded hydropower stations*

Chun-ze ZHANG （張春澤）， Yong-guang CHENG （程永光）， Jia-yang WU （吳家陽）， Wei Diao （刁偉）

State Key Laboratory of Water Resources and Hydropower Engineering Science， Wuhan University， Wuhan 430072， China， E-mail:zhangchunze@whu.edu.cn

?

Lattice Boltzmann simulation of the open channel flow connecting two cascaded hydropower stations*

Chun-ze ZHANG （張春澤）， Yong-guang CHENG （程永光）， Jia-yang WU （吳家陽）， Wei Diao （刁偉）

State Key Laboratory of Water Resources and Hydropower Engineering Science， Wuhan University， Wuhan 430072， China， E-mail:zhangchunze@whu.edu.cn

This paper examines the feasibility and the efficiency of a multiple-relaxation-time lattice Boltzmann model （MRT-LBM）for simulating open channel flows in engineering practice. A MRT-LBM scheme for 2-D shallow water flows taking into account of the bed slope and the friction is proposed. The scheme’s reliability is verified by benchmark problems and the simulation capability is improved by implementing the scheme on a graphic processing unit （GPU）. We use the method to analyze the flow characteristics in the connecting open channel of two cascaded hydropower stations. The flow fields and parameters such as the water depth， the flow rate， and the side-weir discharge， under different operating conditions， are analyzed. The factors affecting the accuracy and the efficiency are discussed. The results are found to be reasonable and may be used as a guidance in the project design. It is shown that the GPU-implemented MRT-LBM on a fine mesh can efficiently simulate two-dimensional shallow water flows in engineering practice.

shallow water equation， lattice Boltzmann method （LBM）， graphic processing unit （GPU）， transient process， cascaded hydropower station

Introduction

For hydropower plants with headrace or tailrace open channels， the transition between different operating conditions generates complicated transient flow problems: the propagation of surge waves along the open channel and the reflection at boundaries， the overflows of over-surged water through side-weirs in the pressure forebay or tailrace ponds， and the flow redistributions among intakes due to different unit operations. These transient flow phenomena should be carefully studied before designing any engineering projects［1，2］. Currently， the most widely adopted approach for simulating open channel flows is to solve the 1-D Saint-Venant equations by the method of characteristics［3］， or the finite difference method［4］， or the finite element method［5］. However， when the multidimensional flow characteristics are prevailing， such as in the side-weir overflow and in the interferential redistribution of flow patterns as discussed in this paper， the use of 2-D or 3-D equations and methods is imperative. A great progress is made in recent years in using 2-D or 3-D methods for open channel flows［6，7］，but the transient processes in hydropower plants remain rarely touched issues［8，9］. This is probably due to the fact that such problems requires very fine meshes and very small time steps， and it is very timeconsuming obtaining an initial steady operating condition before simulating transient processes.

The lattice Boltzmann method （LBM） is an alternative numerical tool for simulating fluid flows. Derived from statistical physics the LBM can simulate flows that are governed by the shallow water equations（SWE） by using the linearized Boltzmann equation. The LBM enjoys many significant advantages such as easy parallelization and convenient implementation of boundary conditions for complex boundaries. At present， the LBM for shallow water equations （LBMSWE） has already been applied to solve many problems: the atmospheric circulation of the northernHemisphere［10］， the wind-driven ocean circulation［11］，the tsunami surge waves［12］， and the flow through Strait of Gibraltar［13］. It has also been applied to solve open channel flows including the 2-D dam break［14］and the surge waves in channels of hydropower plants［8］. However， few attempts were made using this method for simulating transient processes in hydropower plants.

This paper simulates the 2-D shallow water flows in a channel that connects two cascaded hydropower stations using the multiple-relaxation-time LBM（MRT-LBM）. The MRT-LBM models the bed slope and the friction by external forcing terms， which is first verified by benchmark tests. To increase the computing ability， the model is then parallelized on a graphic processing unit （GPU） through the Compute Unified Device Architecture （CUDA） platform［15］. Finally， simulations on sufficiently fine grids are performed for the open channel flow connecting two cascaded hydropower plants and careful investigations are made to deal with the interesting transient flow phenomena.

1. Lbm for shallow water flows

1.1 Shallow water equations

The 2-D shallow water equations （SWE） can be expressed as［16］

The four terms on the right hand side of Eq.（2） are the influence of the bed gradient， the bed shear stress， respectively. In Eq.（2），is the bed elevation，is the water density，is the bed shear stress indirection， which is defined by the Manning coefficientand the depth-averaged （Manning’s coefficient at the bed） velocities

1.2 Lattice Boltzmann model

This paper uses the LBM to numerically solve the flows governed by the SWE. The lattice Boltzmann equation may be expressed as［17］

The local equilibrium distribution functions for the LBM method corresponding to the SWE can be chosen as

The multi-relaxation time （MRT）［16］collision operator is applied. Compared with the single-relaxation time LBM （SRT-LBM）， the MRT-LBM incorporates more fluid physical information and is more stable and accurate.

The particle distribution functions are transformed into moment spaces before the relaxation step:， wheredenotes the transformation matrix from the distribution functions to the momentsis the equilibrium momentswhich is composed of

The diagonal collision matrixcontains the relaxation parameters， which are partly related to the kinematic viscosityνvia

The remaining relaxation factors can be selected in the range. The optimal values for these parameters depend on the specific geometry， the initial and boundary conditions of the system. In our study， 1.0 is used for the parameters except

1.3 Body force terms

In the SWE， the influences of the bed gradient and the bed shear stresses are modeled as the body forces （Eq.（2））. These external forcing terms may be added into the equations of the particle distribution functions in the LBM. We choose the second-order method for introducing the forcing terms［18］， in which the values of the forces at the mid-point between nodeand nodeare used

The forcing term added to the equation of the particle distribution functions （i.e.， Eq.（3）） is

The term of the bed gradient is

which may be simplified to

1.4 Boundary conditions

Boundary conditions play a significant role in ensuring the accuracy and the stability of the LBM. There are two types of boundary conditions used in our study as described below:

（1） Periodic boundary condition

This boundary condition is used to minimize the influence of-direction parameters for the entire system in calculating a 1-D problem with a 2-D model（Examples 2.1 and 2.2）.

（2） Inflow and outflow boundary conditions

Therefore， the unknown distributions can be defined as

1.5 Stability requirements

Like any other numerical methods， the LBM for 2-D shallow water flows must satisfy the stability requirements， which can be expressed as［18］

2. Model verification

Two benchmark problems are simulated to verify the reliability and the performance of the parallel LBM algorithm implemented on the GPU hardware. One test involves a steady flow while the other involves an unsteady tidal wave flow. The accuracy of the simulations is demonstrated by comparing the computational results with the analytical solutions. The performance of the parallel MRT-LBM and SRT-LBM algorithms are compared.

Fig.1 Velocities and free surface profiles for the steady flow over a bump

2.1 Steady flow over a hump

This case is used to investigate the LBM’s ability of recovering the correct steady solution for shallow water flows on a non-flat boundary. In Fig.1， a bump in the flow over a 25 m long channel is defined by

The initial conditions are given byand. On the channel boundaries， the dischargeis used as an inlet boundary condition and the water depthis used as an outlet boundary condition. Figure 1 shows the bed， the water surface， and the velocityprofiles. Above the bump， the depthaveraged velocity reaches its maximum at the peak of the convex profile， and the water surface reaches its minimum at the nadir of the concave profile. The analytical solutions can be obtained from the mass conservation and the Bernoulli’s equations. The analytical values at the bump topare the minimum water leveland the maximum depth-averaged velocity

We first simulate this 1-D problem with a 2-D LBM model. Following references［13，18］， the values for the macroscopic reference velocity and the relaxationtime are set asand， respectively. To analyze the influence of the mesh resolution， several meshes of 512 to 4096 nodes in the longitudinal direction （corresponding to a grid spacing in the range fromto） are selected with the number of nodes along the transversal direction being fixed at 32. The corresponding time steps are in the range fromto. The macroscopic reference velocity is kept constant for all mesh sizes. Both the single precision and the double precision are used for all variables in the GPU implementation to compare the efficiency of the implementation.

Figure 2 shows the results of the water level over the bump using different meshes and precisions. It is clear that for the single precision case， the results on the fine mesh are not closer to the analytical solution than the results on the coarse mesh.

Fig.2 Water level profiles above the bump by MRT-LBM for comparing the influences of GPU digital precision and mesh resolution

A similar finding was also reported in Ref.［12］. The reason for this is that the computations on a finer mesh result in a greater truncation error accumulations for the single precision case. However， for the double precision case， all meshes considered have a good convergence， and it is evident that the finer meshes produce results closer to the analytical solution. Using the double precision is more accurate but requires more computing resources. However， in view of the fact that the relative errors for single precision variables are all smaller than 10-4and the requirements of most engineering applications can be met， we use the single precision in the GPU programming in the remainder of the study.

Table 1 Performance （in MNUPS） of SRT-LBM and MRTLBM in the simulations for different meshesNX× NY

To compare the difference of performances between the SRT-LBM and the MRT-LBM， they both are implemented on a single Tesla C1060 chip using the single precision with simulations on different meshes. Table 1 lists the values of the performance parameter MNUPS （million node updates per second）. It can be seen that the two LB models are almost equivalent in terms of performance on the GPU chip. The volume of data transmission between the host device and the kernel functions for the two LB models are nearly the same， therefore， the bandwidths are also the same. Compared with the SRT-LBM， the MRTLBM requires just a few additional fused multiply-add operations on the GPU chip. If the same algorithms are implemented on a CPU chip， the MRT-LBM would need an extra 15% computational cost as compared to the SRT-LBM［16］. The maximum performance for our case is lower than that previously reported of 500 to 600 MNUPS in reference［14］. This is because the additional non-local memory access and the computations are needed in calculating the bottom gradient forcing term in our case.

2.2 Unsteady tidal wave flow

A tidal wave flow is simulated in a frictionlesschannel of length. The bed elevation （see Fig.3） is given by

The initial conditions for the water depth and velocity are

Fig.3 Numerical and analytical free surfaces for the tidal wave flow at time

Fig.4 Velocity distribution profiles for the tidal wave flow on different meshes at time

Fig.5 Histories of errorand relative error6.84 m）

Fig.6 Schematic diagram of the connective structure of the cascade hydropower station

At the channel inlet and outlet， the water depth and velocity are prescribed as follows

The analytical solutions of this problem may be expressed as［20］

This 1-D problem is calculated by using a 2-D model on a GPU chip. The lattice spacing is refined fromto. Because the macroscopic reference velocityis kept constant for all meshes， the corresponding time step size ranges fromdown to. The relaxation time for the MRTLBM is

Figure 3 shows the simulated free surface on a mesh of， which is consistent with the analytical solution for a very long time， wheredenotes the bed elevation. Figure 4 presents the velocity profiles along the flow direction at the time9 220 s. The results obtained from finer meshes ofandshow excellent agreement with the analytical solution. However， the errors on a coarser mesh ofare relatively large. We define the maximum errorand the maximum relative error， in whichis the LBM result andis the asymptotic analytical solution.

Figure 5 shows the histories of the errorsandfor the mesh. The curves fluctuate but see no tendency to increase the extent of fluctuations with time. The periodical fluctuations are due to the periodicity of the water height specified at the inflow boundary. In the whole flow field， the maximum erroris smaller than 0.0015 m/s， and the relative erroris smaller than 0.5%. These two observations indicate that the accuracy of the single precision GPU computations by the LBM is sufficient for the simulation of unsteady flows in engineering practice.

With much less requirement for memory， the single precision GPU computing can achieve approximately eight times of the performance of the double precision equivalent. Therefore， we suggest the use of the single precision GPU computing for simulations of practical engineering problems.

3. Flows in the connecting open channel between two cascated hydropower stations

3.1 Description of the problem

We simulate flows in the connecting open channel between the tailrace pool of the hydropower station 1 （HS-1） and the pressure forebay of the hydropower station 2 （HS-2）， and study the propagation and reflection of open channel surges under different operating conditions. The functions of the tailrace pool，the connecting channel， and the pressure forebay are to collect the tailrace water from the HS-1’s draft tubes，delivering the flow smoothly to the HS-2’s intakes，and distributing the water uniformly for different intakes， respectively. The pressure forebay is normally a gradually broadened and deepened channel. The forebay’s front end has several intakes collocated for supplying the water through penstocks to turbines with a side overflow weir generally placed in the side bank for draining the extra water. The large body of the water in the forebay can reflect the water hammer pressure from the penstocks and suppress the severe water-level fluctuation in the channel. During and after the load release and the load increase， there are obvious surge wave propagations and reflections along the connecting channel with the water overflowing through the side weir.

Table 2 The typical cases of the cascaded hydropower stations to be simulated

The vertical view and the longitudinal sectionview of the connecting channel are shown in Fig.6. The channel is in the shape of a rectangle of 440 m in total length and 30 m in width along its entire length. The channel is composed of three sections: the first section includes the tailrace pool of the HS-1 and the inclined channel of 240 m in length with the bed slope 1:4000， the second section is the horizontal channel of 100 m in length with the bed slope 0， the final section is the pressure forebay of the HS-2 with a length of 100 m and a bed slope of 9.16:100. The four draft tube outlets of the HS-1 and the four penstock intakes of the HS-2 are located in the upstream end of the tailrace pool and the downstream end of the forebay，respectively. The width for all draft tube outlets and penstock intakes is 5 m. In the model， the relative elevation is used and the intake floor of the HS-2 is taken as the datum point （0 m in elevation）. Accordingly， the elevation of the draft tube of the HS-1 is 9.16 m. A side weir of 50 m in length and 14.33 m in crest elevation is set on the left bank of the pressure forebay of the HS-2. The minimal， normal， and maximal levels of the forebay are 12.67 m， 14.13 m and 15.86 m， respectively. The bed roughness of the whole channel is. The rated discharge of each turbine is 32.2 m3/s.

3.2 Simulation conditions and treatments

The rationality of the design of the cascaded hydropower stations is analyzed through simulations of the open channel waves under the three typical operating conditions specified in Table 2. The selection of these operating conditions is based on the following considerations. For Case 1， since four turbines of the HS-2 release all loads， the tailrace water of the HS-1 will overflow from the side weir； therefore，the maximum discharge and the history of the water levels at probe points should be analyzed. For Case 2，the water level of the entire channel will decline gradually after four turbines of HS-1 release all loads，therefore， the time when the HS-1 releases all loads and the time when the water level of the pressure forebay reaches the allowable minimum water level will be determined. For Case 3， since the units #3 and #4 take on the load in succession， the water level in the pressure forebay will continue to decline over a certain period； therefore， accompanied by a certain extent of the water level fluctuation， the decline of the water level will be monitored and the submerged depth will be analyzed. The MRT-LBM is adopted to simulate all cases， which are implemented on a single Tesla C1060 chip using the single precision. The simulation processes include the following two steps.

First， the initial steady operating conditions are simulated. For a given operating condition， we specify the velocity boundary conditions on the draft tube outlets and the penstock intakes， then conduct the simulation until the fluctuations of all probed parameters are within the allowable tolerance and a steady flow state is attained. The allowable tolerance of the probed parameter is set as， whereis the change of the probed parameter in every 100 steps.

Second， the transient processes are simulated. The simulation starts from the initial flow by imposing the velocity profiles on the draft tube outlets and the penstock intakes and it continues until the detailed transient histories of the water level and the velocity are obtained.

Because of the narrow and long shape of the channel， we specify a slip boundary on both banks of the channel to reduce the dissipation from the boundaries. The velocities specified on the boundaries are calculated based on the prescribed discharge and the water level. For the side weir， when the water level is higher than the crest of the weir， the overflow condition will be applied， otherwise， the slip boundary willbe applied. As the overflow process of the side weir cannot be simulated directly with a 2-D LBM， we treat the boundary as a velocity Dirichlet boundary， by imposing the velocity using the empirical weir flow formula

To monitor the histories of the parameters， we assign 15 probe points （Fig.6） in the upstream end of the HS-1’s tailrace pool， the middle reach of the channel， and the rear end of the HS-2’s pressure forebay， and these probe points are denoted as U1 to U5，M1 to M5， and D1 to D5， respectively.

3.3 Results

This section evaluates the propagation and reflection of the open channel waves after applying changes in the operating conditions.

Fig.7 Histories of the parameters for different meshes in Case 1

3.3.1 Lattice dependence of solutions

To examine the lattice dependence of solutions，we run the LBM code on three different meshes with lattice size0.33 m， 0.2 m and 0.125 m to simulate Case 1（shown in Table 2）. Figure 7 presents the histories of the side weir discharge and the water level of the probe points of U3， M3， and D3 as a consequence of sudden load release of the HS-2， wheredenotes the discharge.

All four units of the HS-1 normally operate after all units of the HS-2 have released the loads， therefore，when the channel waves are damped and approach a steady state， the discharge of the side weir should be equal to 128.8 m3/s， which is the sum of the discharges of the HS-1’s four turbines. It is evident that the results obtained using the two fine meshes of0.2 m and 0.125 m are reasonable and in good agreement with one another， but the results obtained with the coarse mesh ofhave relatively large deviations. The stable discharges derived from the two fine mesh simulations are 126.55 m3/s and 128.02 m3/s with relative errors of 1.74% and 0.61%， respectively. As for the histories of the water level at the probe points U3， M3， and D3， the maximum difference between the two fine meshes is 0.026 m， with a relative difference of about 0.19%. Both fine mesh simulations have good accuracy. In view of the fact that the computing cost with the meshis just 24% more than that with the mesh， we applymesh for all following simulations， and maintain the reference velocityand the simulation time step

Figure 8 shows the longitudinal water surface profiles for the cases of initial operating conditions simulated， wheredenotes the surface profile. The shadow area is the bed of the channel. Due to the negligible degree of variation of the bottom slope， the water surface is almost horizontal along the first 340 m of the channel. The bed slope of the pressure forebay is classified as a steep slope （the gradient is 9.16:100 and the critical depth）， therefore， the flow is diminished longitudinally and the surface rises slightly.

Fig.8 Longitudinal water surface profiles for the cases of initial conditions

3.3.2 Case 1-determination of the maximum overflow discharge

During the transient processes in Case 1， a surge wave is first generated in the pressure forebay due to the decrease of the turbine discharge， from 100% to 0% after the load release of the HS-2’s units atThis wave propagates upstream along the channel and the water level continues to rise. When the water level in the pressure forebay rises over the side weir crest，overflow begins. After the overflow， the level gradually approaches a stable state. The counterflow wave propagates continuously until the wavefront arrives at the draft tube outlets around

At that moment， the wave is reflected by the upstream walls， and the superposition of the reflected wave and the coming wave causes the water level at the probe point U3 to attain its first peak value. Then the reflected wave begins propagating downstream，which leads to a further rising of the water level along the channel. When the reflected cocurrent flow wave appears at the inlet of the pressure forebay again， the water level at the probe point D3 and the discharge of the side weir attain their first peak values. Subsequently， the waves continue to propagate in the channel and are reflected at the upstream and downstream boundaries， and the discharge of the side weir and the water levels drastically change over every cycle （about 120 s） until after 600 s. During the process， the final stable water level is 15.43 m， which is lower than the highest design-specified level of the pressure forebay. Despite of the load release by all HS-2 units， the excessive water overflow from the side weir ensures that the tailrace water level of the HS-1 maintains at 15.37 m， which can mitigate the adverse impact of the head fluctuation on the output of the power of turbines. It is evident that the side weir is reliable and stable under extremely adverse operating conditions.

Fig.9 Histories of water levels at different probe points in Case 2

3.3.3 Case 2-determination of the safe reaction time

In general， if the water level of the pressure forebay of the HS-2 is lower than the minimum safety level， air-entraining vortices will occur， which will pose a threat to the safety of the hydropower station. Therefore， the units of the HS-2 should be shut down before this occurs. The time interval between the time when the HS-1 releases the loads and the time when the HS-2 is shut down safely is called the security reaction time interval. The Case 2 （shown in Table 2）is analyzed to determine this time interval for the HS-2. Figure 9 shows the water level histories at the probe points in Case 2. As the units of the HS-1 are shut down within 10 s， the water level at the channel upstream end declines linearly and the falling wave propagates downstream along the channel. Figure 10 shows the histories of the wave front propagation and reflection. The wave front arrives at the rear wall of the pressure forebay after approximately 60 s， and subsequently the forebay water level begins to decline linearly. After 187 s when the water wave has propagated approximately 1.5 cycles in the channel， the water level of the pressure forebay reaches the allowable minimum level. Therefore， if the HS-1 fully releases the load， a safe operation can be ensured by shutting down all HS-2 units within 180 s.

Fig.10 Histories of waterfront propagation and reflection in Case 2

Fig.11 Histories of water levels at probe points and the discharge of the side weir in Case 3

3.3.4 Case 3-determination of the intensity of water level fluctuation

For units running in a hydraulically connectedsystem， if the working condition is changed for some units for whatever reason， the generated fluctuation in the water level and the discharge will affect the working parameters of the non-action units. This phenomenon is called the hydraulic interference. The Case 3（shown in Table 2） is simulated to analyze the interference of the successive load acceptance of the units #3 and #4 on the operation of the units #1 and #2. Initially， all four HS-1 units and the units #1 and #2 of the HS-2 operate normally， and the side weir overflow discharge is almost equal to the discharge of the two turbines. After units #3 and #4 take on the load in succession， the water level in the pressure forebay declines， the flow ceases coming from the side weir，and the water level gradually approaches a steady state. Figure 11 shows the water level histories at the probe points and the discharge of the side weir in Case 3. When the #3 turbine starts to take on the load， the forebay water level declines linearly. After 10 s， the discharge through the #3 intake attains 100%， meanwhile， the increasing trend of the discharge terminates instantly， and this leads to a slight water level fluctuation. Additionally， because the #4 turbine also begins to take on the load， the water level in the pressure forebay continues to decrease. When the discharge of the #4 intake stops increasing， at approximately 20 s，the water level in front of the inlet fluctuates slightly again. Figure 12 shows the water level at the pressure forebay at different moments of the simulation （the same color scales are used）. The wave continues to propagate in the channel， and at approximately 110 s，the water level near the side weir front is lower than the weir crest elevation， so that the flow ceases coming from the weir， which generates a counterflow wave propagating upstream due to the inertia of the water. At approximately 130 s， the side weir stops working and the water level of the channel ceases to decrease. The counterflow wave continues to propagate in the channel but the wave amplitude gradually attenuates.

In the course of the first cycle （approximately 120 s）， the overall trend of the water level variation is similar to Case 2. Some slight fluctuation generated in the transient process has a trivial effect on the water level in front of the #1 and #2 intakes in the pressure forebay. Over the entire duration of the transient process， the water level in the pressure forebay is maintained in a safe operating range and no severe upheaval or incessant fluctuation occur. It is evident that，the steady operation of units #1 and #2 can be guaranteed in Case 3.

Fig.12 Contours of water level at pressure forebay in different moments

4. Conclusions

The feasibility and the efficiency of the MRTLBM are tested in simulating open channel flows in engineering practice and its application in the transient processes of two cascaded hydropower stations. After assessing the accuracy and the efficiency of the method by two benchmark problems， we apply this approach to analyze the practical engineering problem successfully. The results show that the histories of all relevant quantities of the flow field can be well captured by using a fine mesh and a small time step and themethod is reasonable and efficient to solve engineering problems in practice. Some conclusions can be drawn as follows:

A single precision computation is efficient and can be used to obtain reasonably good results for engineering problems such as the open channel flow in our work. Although the results of the single precision are not monotonously convergent to the analyticcal solutions as those of the double precision when the mesh refinement is performed， the relative errors of most quantities are smaller than 10-4， which meets the accuracy requirement of most hydraulic engineering applications.

The performance of the MRT-LBM is similar to that of the SRT-LBM if it is implemented by GPUs on fine meshes. We suggest to use the MRT-LBM in simulations for hydraulic problems in practice because it is more stable. The results obtained for the wave propagation and reflections in the connecting channel show that the simulations using a fine mesh and a small time step can capture many complex and detailed fluctuation features， which are informative and useful for an in-depth analysis of underlying physical mechanisms.

Some improvement and extension of the method，for example， the treatment of the trapezoidal or compound channel section， can be made in the future to increase the capability of the method.

References

［1］ YANG Jian-dong， ZHAO Kun and LI Lin et al. Analysis on the causes of units 7 and 9 accidents at Sayano-Shushenskaya hydropower station［J］. Journal of Hydroelectric Engineering， 2011， 30（4）: 226-234（in Chinese）.

［2］ ZHANG Xiao-xi， CHENG Yong-guang. Simulation of hydraulic transients in hydropower systems using the 1-D-3-D coupling approach［J］. Journal of Hydrodynamics， 2012， 24（4）: 595-604.

［3］ FUSMBA M. Contribution on transient flow modelling in storm sewers［J］. Journal of Hydraulic Research， 2002，40（6）: 685-693.

［4］ VUGDELIJA M.， STPKANOVI? ?. Determination of a time step interval in hydraulic systems transients simulation［J］. Advances in Engineering Software， 2000， 31（2）: 143-148.

［5］ CHEN J.， STEFFLER P. and HICKS F. Conservative formulation for natural open channels and finite-element implementation［J］. Journal of Hydraulic Engineering，ASCE， 2007， 133（9）: 1064-1073.

［6］ HAN S.， BIRON P. M. and RAMAMRTHY A. S. Threedimensional modelling of flow in sharp open-channel bends with vanes［J］. Journal of Hydraulic Research，2011， 49（1）: 64-72.

［7］ ZENG J.， CONSTANTINESCU G. and WEBER L. 3D calculations of equilibrium conditions in loose-bed open channels with significant suspended sediment load［J］. Journal of Hydraulic Engineering， ASCE， 2010， 136（9）: 557-571.

［8］ ZHANG Ming-liang， XU Yuan-yuan and QIAO Yang et al. Numerical simulation of flow and bed morphology in the case of dam break floods with vegetation effect［J］. Journal of Hydrodynamics， 2016， 28（1）: 23-32.

［9］ CHENG Yong-guang， ZHANG Hui. Numerical simulations of 2-D hydraulic transients using lattice Boltzmann method［J］. Journal of Hydraulic Engineering， 2001， （2）: 32-37（in Chinese）.

［10］ FENG Shi-de， ZHAO Ying and TSUTAHARA M. et al. Lattice Boltzmann model in rotational flow field［J］. Chinese Journal of Geophysics， 2002， 45（2）: 164-170（in Chinese）.

［11］ ZHONG L.， FENG S. and GAO S. Wind-driven ocean circulation in shallow water lattice Boltzmann model［J］. Advances in Atmospheric Sciences， 2005， 22（3）: 349-358.

［12］ JANSSEN C. F.， GRILLI S. T. and KRAFCZYK M. Efficient simulations of long wave propagation and runup using a LBM approach on GPGPU hardware［C］. Proceedings of the 22nd （2012） Offshore and Polar Engineering Conferences. Rodos， Greece， 2012.

［13］ TH?MMES G.， SEA?D M. and BANDA M. K. Lattice Boltzmann methods for shallow water flow applications［J］. International Journal for Numerical Methods in Fluids，2007， 55（7）: 673-692.

［14］ ZHANG Chun-ze， CHENG Yong-guang and LI Yongchang. Lattice Boltzmann model for two-dimensional shallow water flows using parallel computing on GPU［J］. Chinese Journal of Hydrodynamics， 2011， 26（2）: 194-200（in Chinese）.

［15］ T?LKE J. Implementation of a lattice Boltzmann kernel using the compute unified device architecture developed by nVIDIA［J］. Computing and Visualization in Science，2010， 13（1）: 29-39.

［16］ HUMIèRESHD.Multiple-relaxation-time lattice Boltzmann models in three dimensions［J］. Philosophical Transactions of the Royal Society of London. Series A: Mathematical， Physical and Engineering Science， 2002，360（1792）: 437-451.

［17］ ROTHMAN D. H.， ZALESKI S. Lattice-gas cellular automata: Simple models of complex hydrodynamics［M］. Cambridge， UK: Cambridge University Press， 2004，5.

［18］ ZHOU J. G. A lattice Boltzmann model for the shallow water equations［J］. Computer Methods in Applied Mechanics and Engineering， 2002， 191（32）: 3527-3539.

［19］ ZOU Q.， HE X. On pressure and velocity boundary conditions for the lattice Boltzmann BGK model［J］. Physics of Fluids， 1997， 9（6）: 1591-1599.

［20］ BERMUDEZ A.， VAZQUEZ M. Upwind methods for hyperbolic conservation laws with source terms［J］. Computers and Fluids， 1994， 23（8）: 1049-1071.

August 22， 2014， Revised December 8， 2014）

* Project support by the National Natural Science Foundation of China （Grant Nos. 11172219， 51579187）， the Specialized Research Fund for the Doctoral Program of Higher Education of China （Grant No. 20130141110013）.

Biography: Chun-ze ZHANG （1986-）， Male， Ph. D.

Yong-guang CHENG，

E-mail:ygcheng@whu.edu.cn