- Email: [email protected]
A lattice Boltzmann model for the shallow water equations
Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539 www.elsevier.com/locate/cma
A lattice Boltzmann model for the shallow water equations J.G. Zhou Department of Computing and Mathematics, Manchester Metropolitan University, Manchester, UK M1 5GD Received 19 December 2000; received in revised form 9 November 2001; accepted 21 February 2002
Abstract A lattice Boltzmann model for the shallow water equations (LABSWE) is developed. It solves the equations with the source terms such as bed slope, bed friction. The model has been verified by solving both steady and unsteady flow problems. Excellent agreement is obtained between numerical predictions and analytical solutions. The results show that the LABSWE is a very competitive method for computational hydraulics or hydrodynamics in terms of computational efficiency and accuracy. 2002 Published by Elsevier Science B.V. Keywords: Lattice Boltzmann method; Shallow water equations; Source terms
1. Introduction The lattice Boltzmann model is evolved from the lattice gas model in order to overcome the difficulties with the lattice gas model [1,2]. It is based on the statistical physics and describes the microscopic picture of particles movement in an extremely simplified way, but on the macroscopic level it gives a correct average description of a fluid. The averaged particle velocities behave in time and space just as the flow velocities in a physical fluid. This provides an indirect way to solve flow equations. The main advantages are that: (1) it consists of simple arithmetic calculations, hence it is easy to program; (2) only one single variable, the microscopic distribution function, is unknown and to be determined; (3) the current value of the distribution function depends only on the previous conditions which is ideal for parallel computations; (4) it is suitable for flows in complex geometry such as flows through porous media because of easy implementation of boundary conditions; and (5) it is easy to simulate complex flows, for example, multiphase flows and flows with variations of boundaries. It is these features that make the lattice Boltzmann method a very promising computational method in different areas [2]. However, most of the researches reported in the literature are limited to the lattice Boltzmann methods for Navier–Stokes equations. On the other hand, the shallow water equations have wide range of applications in ocean, environmental and hydraulic engineering, for instance, tidal flows in estuary and coastal regions, river, reservoir and open channel flows. There are many computational methods for the shallow water equations such as finite difference method, finite volume method, finite element method and Godunov-type method [3–5]. Moreover, it
E-mail address: [email protected] (J.G. Zhou). 0045-7825/02/$ - see front matter 2002 Published by Elsevier Science B.V. PII: S 0 0 4 5 - 7 8 2 5 ( 0 2 ) 0 0 2 9 1 - 8
3528
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
may be noted that for a Godunov-type method, the inclusion of source terms, e.g. bed slope term, often causes difficulties in solution of the equations [6]; while, for a non-Godunov type method, it is straightforward to add source terms in the equations. All of these methods are developed based on direct solutions of the shallow water equations. Therefore, in this paper, we develop a lattice Boltzmann model for the shallow water equations (LABSWE). It provides an alternative method for the solution of the shallow water equations with or without source terms. The model has been verified and its accuracy has been demonstrated by solving both steady and unsteady flow problems. 2. Lattice Boltzmann model In this section, first of all, the 2D shallow water equations are introduced. Then the lattice Boltzmann equation is constructed for the solution of the shallow water equations. Finally, the boundary conditions are described. 2.1. Shallow water equations The 2D shallow water equations with source terms of bed slope and bed friction may be written in tensor notation as [3] oh oðhui Þ þ ¼ 0; ð1Þ ot oxi oðhui Þ oðhui uj Þ o h2 ozb sbi o oui þ ¼ g þ hm ; ð2Þ gh ot oxj oxi 2 oxi q oxj oxj where i and j are indices and the Einstein summation convention is used, i.e. repeated indices mean a summation over the space coordinates; xi is the Cartesian coordinate; h the water depth; t the time; ui the depth-averaged velocity component in i direction; zb the bed elevation above datum; g ¼ 9:81 m/s2 the gravitational acceleration; q the water density; m the kinematic viscosity; sbi the bed shear stress in i direction defined by the depth-averaged velocities pffiffiffiffiffiffiffiffi sbi ¼ qCb ui uj uj ; ð3Þ in which Cb is the bed friction coefficient, which may be either constant or estimated from Cb ¼ g=Cz2 , where Cz ¼ h1=6 =nb the Chezy constant (nb the Manning’s coefficient at the bed). 2.2. Lattice Boltzmann equation According to the theory of the lattice Boltzmann method, it consists of two steps: a streaming step and a collision step. In the streaming step, the particles move to the neighbouring lattice points which is governed by Dt fa0 ðx þ ea Dt; t þ DtÞ ¼ fa ðx; tÞ þ eai Fi ðx; tÞ; ð4Þ Na e2 where fa is the distribution function of particles; fa0 the value of fa after the streaming; e ¼ Dx=Dt; Dx the lattice size; Dt the time step; Fi the component of the force in i direction; ea the velocity vector of a particle in the a link and Na is a constant, which is decided by the lattice pattern as 1 X Na ¼ 2 eai eai : ð5Þ e a
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3529
In the collision step, the arriving particles at the points interact one another and change their velocity directions according to scattering rules, which is expressed as fa ðx; tÞ ¼ fa0 ðx; tÞ Xa ½fa0 ðx; tÞ faeq ðx; tÞ ;
ð6Þ
faeq
where is the local equilibrium distribution function and Xa the collision operator which controls the speed of change in fa from collision. Usually, with the BGK approximation [7] under which the collision operator Xa is replaced with the single relaxation time s, these two steps can be combined into the following lattice Boltzmann equation: 1 Dt fa ðx þ ea Dt; t þ DtÞ fa ðx; tÞ ¼ ðfa faeq Þ þ eai Fi ðx; tÞ; s Na e2
ð7Þ
which is the common form of a lattice Boltzmann model. The physical variables, the water depth h and the velocity ui , are defined in terms of the distribution function as X 1X h¼ fa ; ui ¼ eai fa : ð8Þ h a a Although there are two main lattice patterns, square lattice and hexagonal lattice, suggested in the literature, the 9-speed square lattice shown in Fig. 1 is frequently used. It can provide an easy way to implement boundary conditions and usually gives more accurate results [8]. Moreover, only using a square lattice, can the force term associated with the bed slope be accurately represented in the model, i.e. the bed elevation zb is simply defined at the centre of a lattice which minimizes error in calculation of its gradient. Hence the 9-speed square lattice is used in the present study (see Appendix A for a lattice Boltzmann model using hexagonal lattice). The velocity vector of particles is defined by 8 0; a ¼ 0; > > > pða 1Þ pða 1Þ < Þ þ~ |e sinð Þ; a ¼ 1; 3; 5; 7; ð9Þ ea ¼ ~ıe cosð 4 4 > p ffiffi ffi p ffiffi ffi > pða 1Þ pða 1Þ > :~ı 2e cosð Þ þ~ | 2e sinð Þ; a ¼ 2; 4; 6; 8: 4 4 With the square lattice it can be shown Na ¼ 6. In order to apply the lattice Boltzmann equation (7) for the solution of the 2D shallow water equations (1) and (2), a suitable local equilibrium function faeq should be decided. In the lattice Boltzmann dynamics, faeq is often assumed to be expressed as a power series in the bulk velocity, i.e. faeq ¼ A þ Beai ui þ Ceai eaj ui uj þ Dui uj dij ;
Fig. 1. 9-Speed square lattice.
ð10Þ
3530
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
where the Kronecker delta function dij is 0; if i 6¼ j; dij ¼ 1; if i ¼ j: The local equilibrium function (10) must satisfy the following conditions: X faeq ¼ h;
ð11Þ
ð12Þ
a
X
eai faeq ¼ hui ;
ð13Þ
eai eaj faeq ¼ 12gh2 dij þ hui uj ;
ð14Þ
a
X a
so that the solution of the lattice Boltzmann equation (7) approaches the solution of the 2D shallow water equations (1) and (2). With the conditions (12)–(14) the coefficients A, B, C and D are determined in a similar way to that of the lattice Boltzmann equation for the Navier–Stokes equations which is detailed by Rothman and Zaleski [9], hence it is not repeated here. After the coefficients are decided, the local equilibrium function reads: 8 5gh2 2h > > 2 ui uj dij ; a ¼ 0; h > 2 > 3e > < 2 6e gh h h h ð15Þ faeq ¼ þ 2 eai ui þ 4 eai eaj ui uj 2 ui uj dij ; a ¼ 1; 3; 5; 7; 2 > 3e 2e 6e 6e > > 2 > > : gh þ h eai ui þ h eai eaj ui uj h ui uj dij ; a ¼ 2; 4; 6; 8: 8e4 24e2 24e2 12e2 By applying the Chapman–Engkog procedure, it can be shown that the solution of the lattice Boltzmann equation (7) with the equilibrium function (15) results in the solution of the shallow water equations (1) and (2) with Fi ðx; tÞ ¼ gh
ozb sbi ; oxi q
m¼
e2 Dt ð2s 1Þ: 6
ð16Þ
It must be pointed that more source terms such as wind shear stress and Coriolis force can easily be added to the lattice Boltzmann equation (7) via Eq. (16). The solution procedure for the LABSWE is now summarized as follows: (1) (2) (3) (4) (5)
given initial water depth and velocity, calculate faeq from Eq. (15), compute fa via the lattice Boltzmann equation (7) with a proper value for the relaxation time s, update the depth and the velocity according to Eq. (8), repeat the steps (2)–(4) until a solution is obtained.
2.3. Boundary conditions There are different boundary conditions in the solution of the shallow water equations. These conditions must be correctly transformed to suitable boundary conditions for the LABSWE. This is briefly described as follows:
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3531
• Inflow, outflow boundary conditions: if the velocity and the depth are known, the unknown distribution function fa at the boundary can be decided with the method described by Zou and He [10]. For example, given velocity at the inflow boundary, after streaming, the unknown f1 , f2 and f8 (see Fig. 1) can be decided as f1 ¼ f5 þ
2hu ; 3e
ð17Þ
f2 ¼
hu f7 f3 þ f6 þ ; 6e 2
ð18Þ
f8 ¼
hu f3 f7 þ f4 þ : 6e 2
ð19Þ
If a zero gradient for velocity or depth at the either inflow or outflow boundary is required, we may set the gradient of the distribution function normal to the boundary to zero. • Solid boundary conditions: we may apply no-slip or slip boundary conditions. For no-slip conditions, the standard bounce-back scheme can be used. For slip conditions, a zero gradient of the distribution function normal to the solid wall can be used.
3. Verification of the model In this section, the LABSWE is used to solve some benchmark problems including both steady and unsteady flow problems. The accuracy is demonstrated by comparing numerical predictions with analytical solutions and available experimental data. 3.1. Steady flow over a bump A 1D steady subcritical flow in a 25 m long channel with a bump defined by 2 0:2 0:05ðx 10Þ ; if 8 < x < 12; zb ðxÞ ¼ 0; otherwise
ð20Þ
is a classical test problem which has been used as a benchmark test case for numerical methods at a workshop on dam-break wave simulations [11]. The problem was also used by Vazquez-Cend on [12] to test their scheme with an upwind discretisation for the bed slope source term. The analytical solution is given by Goutal and Maurel [11]. The global relative error R is defined by vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u uX hn hn1 2 i i R¼t ; ð21Þ n h i i where hn and hn1 are the local water depths at the current and previous time levels. The convergence criterion for a steady solution is defined as R < 5 106 . In the numerical computations, the discharge per unit width of q ¼ 4:42 m2 /s was imposed at the inflow and h ¼ 2 m was specified at the downstream end. e ¼ 15 m/s and s ¼ 1:5. The 1D flow can be guaranteed in the 2D code by specifying slip boundary conditions at side walls. To test the effect of the lattice size on the solutions, three lattices, 125 50, 250 50 and 500 50 corresponding respectively to Dx ¼ Dy ¼ 0:2 m, Dx ¼ Dy ¼ 0:1 m and Dx ¼ Dy ¼ 0:05 m, were used in the initial computations. Steady state solutions were reached after 20,000 steps. Comparison between numerical results based on different lattice sizes
3532
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
are depicted in Fig. 2, showing that there is little difference. Further quantitative comparison indicates that there is small difference of the results between lattice Dx ¼ 0:2 m and Dx ¼ 0:1 m, but the results are almost the same for both Dx ¼ 0:1 m and Dx ¼ 0:05 m. This suggests that the results based on lattice with Dx 6 0:1 m provide lattice-independent solutions; hence the solution obtained with Dx ¼ 0:05 m is presented here. Fig. 3 shows the profile of the water surface along the channel, in which there is a surface drop because the flow is subcritical. Comparison of the result with the analytical solution is plotted in Fig. 4, showing an
Fig. 2. Subcritical flow: effect of the lattice size on solutions.
Fig. 3. Subcritical flow: profile of the water surface.
Fig. 4. Subcritical flow: comparison of the water surface.
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3533
excellent agreement. The quantitative comparison with the analytical solution indicates that the maximum relative error for the water depth is smaller than 0:25%. In order to test the conservative property of the model, the numerical solution of the discharge (mass) is shown in Fig. 5. Also the quantitative comparison with the theoretical discharge shows that the maximum relative error is smaller than 0:15%. This suggests that the model is accurate and conservative. 3.2. Tidal wave flow Tidal waves often occur in coastal engineering. Here we consider the test problem that Bermudez and V azquez [13] used for verification of an upwind discretisation of the bed slope source terms. This is a one dimensional problem with bed topography defined by (see Fig. 6)
40x 4x 1 H ðxÞ ¼ 50:5 10 sin p ; ð22Þ L L 2
Fig. 5. Subcritical flow: comparison of the discharge.
Fig. 6. Tidal flow: comparison of the water surface at t ¼ 9117:5 s.
3534
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
where L ¼ 14,000 m is the channel length and H ðxÞ the partial depth between a fixed reference level and the bed surface, hence zb ðxÞ ¼ H ð0Þ H ðxÞ. The initial and boundary conditions are
and
hðx; 0Þ ¼ H ðxÞ;
ð23Þ
uðx; 0Þ ¼ 0
ð24Þ
hð0; tÞ ¼ H ð0Þ þ 4 4 sin p
4t 1 þ 86;400 2
;
uðL; tÞ ¼ 0:
ð25Þ ð26Þ
Under these conditions, the tidal wave is relatively short and an asymptotic analytical solution is given by Bermudez and V azquez [13] as
hðx; tÞ ¼ H ðxÞ þ 4 4 sin p
4t 1 þ ; 86;400 2
ðx 14;000Þp 4t 1 cos p uðx; tÞ ¼ þ : 5400hðx; tÞ 86;400 2
ð27Þ ð28Þ
To achieve a lattice-independent solution, three lattices of 600, 800 and 1000 were used. e ¼ 200 m/s and s ¼ 0:6. Eqs. (23)–(26) were used as the initial and boundary conditions. Fig. 7 shows the comparison of the results based on the lattices. Difference in the results with 600 and 800 lattices clearly exists, but the results based on 800 and 1000 lattices are almost the same. Hence the results with 800 lattices are further described here. A comparison of the numerical results with the asymptotic analytical solutions at t ¼ 9117:5 s is shown in Figs. 6 and 7, respectively, showing good agreement. The quantitative comparison with the analytical solution is carried out. For the water depths, the maximum relative error is smaller than 1%. For the velocity, the maximum relative error is smaller than 5% when x 6 11,445 m; when x > 11,445 m, the
Fig. 7. Tidal flow: comparison of the velocity uðx; tÞ at t ¼ 9117:5 s.
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3535
relative error may be larger than 5% and it becomes big with x, which is due to the fact that the zero boundary condition was used for velocity at the downstream end, leading to u ! 0 while x ! 14,000 m. Consequently, the velocity becomes very small such that the relative error is no longer an appropriate measurement for accuracy; instead the absolute error is usually used. For this case, the further calculation shows that the absolute error is smaller than 0:006 m/s for all the velocities with the relative error larger than 5%. Since the analytical solution is obtained by use of the asymptotic analysis, such agreement may be considered to be excellent. This confirms that the model is accurate for unsteady shallow water flow problems. The present method can provide solution of the same accuracy as that reported by Bermudez and Vazquez who used a complex upwind discretisation for the bed slope source terms [13]. 3.3. Flow around a cylinder Flow around a cylinder is a classical problem in fluvial hydraulics. In order to test the capability of the present method, a cylinder with radius of 0.11 m situated in the centre of a channel is used, which is the same as Test 1 used by Yulistiyanto et al. in the numerical and experimental investigations [14]. The channel is 4 m long and 2 m wide (see Fig. 9). The discharge is Q ¼ 0:248 m3 /s; outflow depth is ho ¼ 0:185 m; the bed slope is ozb =ox ¼ 6:25 104 in the flow direction; and the Manning’s coefficient is nb ¼ 0:012. In the computation, 600 300 square lattices were used. Dx ¼ Dy ¼ 0:00667 m, Dt ¼ 0:00145 s and s ¼ 1:982. No-slip boundary condition is used along the wall of the cylinder; slip boundary condition is used for the side walls; the depth at the outflow is set to h ¼ ho ; zero gradient for the depth is specified at the inflow boundary; the inflow velocity u is decided based on the discharge; and v ¼ 0 at the inflow boundary. After 40,000 steps, a steady state solution was reached. Comparison of the depths along the centerline of the channel between the computation and the experimental data is shown in Fig. 8. The calculation indicates that the maximum relative error for the depth is 9.2%; hence agreement is reasonably good. The contours of the water depths are plotted in Fig. 9, showing a similar pattern to Fig. 10 obtained in the experiment. Fig. 11 shows the flow characteristics in the vicinity of the cylinder, which is again similar to Fig. 12 given by Yulistiyanto et al. [14] with a different numerical method. In fact, there is a relation between the wake length Lw and the cylinder diameter D, i.e. Lw 1:3D (see Fig. 12 for definitions). This is also in good agreement with that obtained by Yulistiyanto et al. [14]. The vectors of the whole flow field is depicted in Fig. 13. However, the difference exists in the region close to the cylinder, which may be due to the use of both no-slip boundary conditions and uniform square lattices in the computations. The streamlines are plotted in Fig. 14. 3.4. Flow over a sudden-expansion channel In this section, flow over a channel with a symmetric sudden-expansion is considered to demonstrate that the method is able to simulate a recirculation in shallow water flows. The entrance channel is 1 m wide and
Fig. 8. Flow around cylinder: comparison of the depth along channel centerline.
3536
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
Fig. 9. Flow around cylinder: contours of the water depth (computation).
Fig. 10. Flow around cylinder: contours of the water depth (experiment) from Ref. [14].
Fig. 11. Flow around cylinder: velocity vectors (computation).
Fig. 12. Flow around cylinder: velocity vectors from Ref. [14].
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3537
Fig. 13. Flow around cylinder: velocity vectors.
Fig. 14. Flow around cylinder: streamlines.
Fig. 15. Sudden-expansion channel: velocity vectors.
2 m long, and the expanded channel is 3 m wide and 4 m long (see Fig. 15). There is no bed slope in the channel, and the bed friction is neglected. In the computation, 120 60 square lattices were used. Dx ¼ Dy ¼ 0:05 m, Dt ¼ 0:025 s and s ¼ 1. No-slip boundary condition is used along the walls. The depth at the outflow is set to h ¼ 0:16 m; zero gradient for the depth is specified at the inflow boundary; the discharge is Q ¼ 0:032 m3 /s; and v ¼ 0 at the inflow boundary. After 10,000 steps, a steady state solution was reached. The velocity vectors are plotted in Fig. 15, showing well-developed circulating flows on both sides. This confirms that the present method is capable of simulating recirculation occurring in shallow water flows.
3538
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
4. Conclusions An accurate LABSWE is presented in the paper. The model is simple and conservative. It is able to solve steady and unsteady flow problems. The source terms can naturally be added to the lattice Boltzmann equation as a force term without any special treatment. This makes the model have potential capability in solving practical flow problems.
Acknowledgements We would like to thank Professor W.H. Graf (Laboratoire de Recherches Hydrauliques, Switzerland) for his providing the experimental data of the flow around a cylinder.
Appendix A. 7-Speed model A lattice Boltzmann model for the shallow water equations using the 7-speed hexagonal lattice shown in Fig. A.1 can similarly be developed as 1 Dt fa ðx þ ea Dt; t þ DtÞ fa ðx; tÞ ¼ ðfa faeq Þ þ 2 eai Fi ðx; tÞ s 3e
ðA:1Þ
with the velocity vector of particles defined by ea ¼
0; ~ıe cosðpða1Þ Þ þ~ |e sinðpða1Þ Þ; 3 3
a ¼ 0; a¼16
ðA:2Þ
and the local equilibrium function faeq expressed as
faeq
8 gh2 h > > < h 2 þ 2 ui uj dij ; e e ¼ 2 > gh h 2h h > : þ 2 eai ui þ 4 eai eaj ui uj 2 ui uj dij ; 2 3e 3e 2e 6e
a ¼ 0; ðA:3Þ a ¼ 1 6:
Fig. A.1. 7-Speed hexagonal lattice.
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3539
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: theory and applications, Phys. Rep. 222 (1992) 45–197. S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Ann. Rev. Fluid Mech. 30 (1998) 329–364. J.G. Zhou, Velocity–depth coupling in shallow water flows, J. Hydr. Engrg., ASCE 121 (10) (1995) 717–724. P.K. Stansby, J.G. Zhou, Shallow-water flow solver with non-hydrostatic pressure: 2D vertical plane problems, Int. J. Numer. Meth. Fluids 28 (1998) 541–563. E.F. Toro, Riemann problems and the WAF method for solving two-dimensional shallow water equations, Phil. Trans. Roy. Soc. Lond. A 338 (1992) 43–68. R.J. LeVeque, Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wavepropagation algorithm, J. Comput. Phys. 146 (1998) 346–365. P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component system, Phys. Rev. 94 (1954) 511–525. P.A. Skordos, Initial and boundary conditions for the lattice Boltzmann method, Phys. Rev. E 48 (1993) 4823–4842. D.H. Rothman, S. Zaleski, Lattice-Gas Cellular Automata, Cambridge University Press, London, 1997. Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Phys. Fluids 9 (1997) 1591– 1598. N. Goutal, F. Maurel (Eds.), Proceedings of the Second Workshop on Dam-break Wave Simulation, HE-43/97/016/B, Departement Laboratoire National d’Hydraulique, Groupe Hydraulique Fluviale, Electricite de France, France, 1997. M.E. V azquez-Cend on, Improved treatment of source terms in upwind schemes for shallow water equations in channels with irregular geometry, J. Comput. Phys. 148 (1999) 497–526. A. Bermudez, M.E. Vazquez, Upwind methods for hyperbolic conservation laws with source terms, Comput. Fluids 23 (1994) 1049–1071. B. Yulistiyanto, Y. Zech, W.H. Graf, Flow around a cylinder: shallow-water modeling with diffusion–dispersion, J. Hydr. Engrg., ASCE 124 (4) (1998) 419–429.
A lattice Boltzmann model for the shallow water equations J.G. Zhou Department of Computing and Mathematics, Manchester Metropolitan University, Manchester, UK M1 5GD Received 19 December 2000; received in revised form 9 November 2001; accepted 21 February 2002
Abstract A lattice Boltzmann model for the shallow water equations (LABSWE) is developed. It solves the equations with the source terms such as bed slope, bed friction. The model has been verified by solving both steady and unsteady flow problems. Excellent agreement is obtained between numerical predictions and analytical solutions. The results show that the LABSWE is a very competitive method for computational hydraulics or hydrodynamics in terms of computational efficiency and accuracy. 2002 Published by Elsevier Science B.V. Keywords: Lattice Boltzmann method; Shallow water equations; Source terms
1. Introduction The lattice Boltzmann model is evolved from the lattice gas model in order to overcome the difficulties with the lattice gas model [1,2]. It is based on the statistical physics and describes the microscopic picture of particles movement in an extremely simplified way, but on the macroscopic level it gives a correct average description of a fluid. The averaged particle velocities behave in time and space just as the flow velocities in a physical fluid. This provides an indirect way to solve flow equations. The main advantages are that: (1) it consists of simple arithmetic calculations, hence it is easy to program; (2) only one single variable, the microscopic distribution function, is unknown and to be determined; (3) the current value of the distribution function depends only on the previous conditions which is ideal for parallel computations; (4) it is suitable for flows in complex geometry such as flows through porous media because of easy implementation of boundary conditions; and (5) it is easy to simulate complex flows, for example, multiphase flows and flows with variations of boundaries. It is these features that make the lattice Boltzmann method a very promising computational method in different areas [2]. However, most of the researches reported in the literature are limited to the lattice Boltzmann methods for Navier–Stokes equations. On the other hand, the shallow water equations have wide range of applications in ocean, environmental and hydraulic engineering, for instance, tidal flows in estuary and coastal regions, river, reservoir and open channel flows. There are many computational methods for the shallow water equations such as finite difference method, finite volume method, finite element method and Godunov-type method [3–5]. Moreover, it
E-mail address: [email protected] (J.G. Zhou). 0045-7825/02/$ - see front matter 2002 Published by Elsevier Science B.V. PII: S 0 0 4 5 - 7 8 2 5 ( 0 2 ) 0 0 2 9 1 - 8
3528
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
may be noted that for a Godunov-type method, the inclusion of source terms, e.g. bed slope term, often causes difficulties in solution of the equations [6]; while, for a non-Godunov type method, it is straightforward to add source terms in the equations. All of these methods are developed based on direct solutions of the shallow water equations. Therefore, in this paper, we develop a lattice Boltzmann model for the shallow water equations (LABSWE). It provides an alternative method for the solution of the shallow water equations with or without source terms. The model has been verified and its accuracy has been demonstrated by solving both steady and unsteady flow problems. 2. Lattice Boltzmann model In this section, first of all, the 2D shallow water equations are introduced. Then the lattice Boltzmann equation is constructed for the solution of the shallow water equations. Finally, the boundary conditions are described. 2.1. Shallow water equations The 2D shallow water equations with source terms of bed slope and bed friction may be written in tensor notation as [3] oh oðhui Þ þ ¼ 0; ð1Þ ot oxi oðhui Þ oðhui uj Þ o h2 ozb sbi o oui þ ¼ g þ hm ; ð2Þ gh ot oxj oxi 2 oxi q oxj oxj where i and j are indices and the Einstein summation convention is used, i.e. repeated indices mean a summation over the space coordinates; xi is the Cartesian coordinate; h the water depth; t the time; ui the depth-averaged velocity component in i direction; zb the bed elevation above datum; g ¼ 9:81 m/s2 the gravitational acceleration; q the water density; m the kinematic viscosity; sbi the bed shear stress in i direction defined by the depth-averaged velocities pffiffiffiffiffiffiffiffi sbi ¼ qCb ui uj uj ; ð3Þ in which Cb is the bed friction coefficient, which may be either constant or estimated from Cb ¼ g=Cz2 , where Cz ¼ h1=6 =nb the Chezy constant (nb the Manning’s coefficient at the bed). 2.2. Lattice Boltzmann equation According to the theory of the lattice Boltzmann method, it consists of two steps: a streaming step and a collision step. In the streaming step, the particles move to the neighbouring lattice points which is governed by Dt fa0 ðx þ ea Dt; t þ DtÞ ¼ fa ðx; tÞ þ eai Fi ðx; tÞ; ð4Þ Na e2 where fa is the distribution function of particles; fa0 the value of fa after the streaming; e ¼ Dx=Dt; Dx the lattice size; Dt the time step; Fi the component of the force in i direction; ea the velocity vector of a particle in the a link and Na is a constant, which is decided by the lattice pattern as 1 X Na ¼ 2 eai eai : ð5Þ e a
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3529
In the collision step, the arriving particles at the points interact one another and change their velocity directions according to scattering rules, which is expressed as fa ðx; tÞ ¼ fa0 ðx; tÞ Xa ½fa0 ðx; tÞ faeq ðx; tÞ ;
ð6Þ
faeq
where is the local equilibrium distribution function and Xa the collision operator which controls the speed of change in fa from collision. Usually, with the BGK approximation [7] under which the collision operator Xa is replaced with the single relaxation time s, these two steps can be combined into the following lattice Boltzmann equation: 1 Dt fa ðx þ ea Dt; t þ DtÞ fa ðx; tÞ ¼ ðfa faeq Þ þ eai Fi ðx; tÞ; s Na e2
ð7Þ
which is the common form of a lattice Boltzmann model. The physical variables, the water depth h and the velocity ui , are defined in terms of the distribution function as X 1X h¼ fa ; ui ¼ eai fa : ð8Þ h a a Although there are two main lattice patterns, square lattice and hexagonal lattice, suggested in the literature, the 9-speed square lattice shown in Fig. 1 is frequently used. It can provide an easy way to implement boundary conditions and usually gives more accurate results [8]. Moreover, only using a square lattice, can the force term associated with the bed slope be accurately represented in the model, i.e. the bed elevation zb is simply defined at the centre of a lattice which minimizes error in calculation of its gradient. Hence the 9-speed square lattice is used in the present study (see Appendix A for a lattice Boltzmann model using hexagonal lattice). The velocity vector of particles is defined by 8 0; a ¼ 0; > > > pða 1Þ pða 1Þ < Þ þ~ |e sinð Þ; a ¼ 1; 3; 5; 7; ð9Þ ea ¼ ~ıe cosð 4 4 > p ffiffi ffi p ffiffi ffi > pða 1Þ pða 1Þ > :~ı 2e cosð Þ þ~ | 2e sinð Þ; a ¼ 2; 4; 6; 8: 4 4 With the square lattice it can be shown Na ¼ 6. In order to apply the lattice Boltzmann equation (7) for the solution of the 2D shallow water equations (1) and (2), a suitable local equilibrium function faeq should be decided. In the lattice Boltzmann dynamics, faeq is often assumed to be expressed as a power series in the bulk velocity, i.e. faeq ¼ A þ Beai ui þ Ceai eaj ui uj þ Dui uj dij ;
Fig. 1. 9-Speed square lattice.
ð10Þ
3530
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
where the Kronecker delta function dij is 0; if i 6¼ j; dij ¼ 1; if i ¼ j: The local equilibrium function (10) must satisfy the following conditions: X faeq ¼ h;
ð11Þ
ð12Þ
a
X
eai faeq ¼ hui ;
ð13Þ
eai eaj faeq ¼ 12gh2 dij þ hui uj ;
ð14Þ
a
X a
so that the solution of the lattice Boltzmann equation (7) approaches the solution of the 2D shallow water equations (1) and (2). With the conditions (12)–(14) the coefficients A, B, C and D are determined in a similar way to that of the lattice Boltzmann equation for the Navier–Stokes equations which is detailed by Rothman and Zaleski [9], hence it is not repeated here. After the coefficients are decided, the local equilibrium function reads: 8 5gh2 2h > > 2 ui uj dij ; a ¼ 0; h > 2 > 3e > < 2 6e gh h h h ð15Þ faeq ¼ þ 2 eai ui þ 4 eai eaj ui uj 2 ui uj dij ; a ¼ 1; 3; 5; 7; 2 > 3e 2e 6e 6e > > 2 > > : gh þ h eai ui þ h eai eaj ui uj h ui uj dij ; a ¼ 2; 4; 6; 8: 8e4 24e2 24e2 12e2 By applying the Chapman–Engkog procedure, it can be shown that the solution of the lattice Boltzmann equation (7) with the equilibrium function (15) results in the solution of the shallow water equations (1) and (2) with Fi ðx; tÞ ¼ gh
ozb sbi ; oxi q
m¼
e2 Dt ð2s 1Þ: 6
ð16Þ
It must be pointed that more source terms such as wind shear stress and Coriolis force can easily be added to the lattice Boltzmann equation (7) via Eq. (16). The solution procedure for the LABSWE is now summarized as follows: (1) (2) (3) (4) (5)
given initial water depth and velocity, calculate faeq from Eq. (15), compute fa via the lattice Boltzmann equation (7) with a proper value for the relaxation time s, update the depth and the velocity according to Eq. (8), repeat the steps (2)–(4) until a solution is obtained.
2.3. Boundary conditions There are different boundary conditions in the solution of the shallow water equations. These conditions must be correctly transformed to suitable boundary conditions for the LABSWE. This is briefly described as follows:
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3531
• Inflow, outflow boundary conditions: if the velocity and the depth are known, the unknown distribution function fa at the boundary can be decided with the method described by Zou and He [10]. For example, given velocity at the inflow boundary, after streaming, the unknown f1 , f2 and f8 (see Fig. 1) can be decided as f1 ¼ f5 þ
2hu ; 3e
ð17Þ
f2 ¼
hu f7 f3 þ f6 þ ; 6e 2
ð18Þ
f8 ¼
hu f3 f7 þ f4 þ : 6e 2
ð19Þ
If a zero gradient for velocity or depth at the either inflow or outflow boundary is required, we may set the gradient of the distribution function normal to the boundary to zero. • Solid boundary conditions: we may apply no-slip or slip boundary conditions. For no-slip conditions, the standard bounce-back scheme can be used. For slip conditions, a zero gradient of the distribution function normal to the solid wall can be used.
3. Verification of the model In this section, the LABSWE is used to solve some benchmark problems including both steady and unsteady flow problems. The accuracy is demonstrated by comparing numerical predictions with analytical solutions and available experimental data. 3.1. Steady flow over a bump A 1D steady subcritical flow in a 25 m long channel with a bump defined by 2 0:2 0:05ðx 10Þ ; if 8 < x < 12; zb ðxÞ ¼ 0; otherwise
ð20Þ
is a classical test problem which has been used as a benchmark test case for numerical methods at a workshop on dam-break wave simulations [11]. The problem was also used by Vazquez-Cend on [12] to test their scheme with an upwind discretisation for the bed slope source term. The analytical solution is given by Goutal and Maurel [11]. The global relative error R is defined by vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u uX hn hn1 2 i i R¼t ; ð21Þ n h i i where hn and hn1 are the local water depths at the current and previous time levels. The convergence criterion for a steady solution is defined as R < 5 106 . In the numerical computations, the discharge per unit width of q ¼ 4:42 m2 /s was imposed at the inflow and h ¼ 2 m was specified at the downstream end. e ¼ 15 m/s and s ¼ 1:5. The 1D flow can be guaranteed in the 2D code by specifying slip boundary conditions at side walls. To test the effect of the lattice size on the solutions, three lattices, 125 50, 250 50 and 500 50 corresponding respectively to Dx ¼ Dy ¼ 0:2 m, Dx ¼ Dy ¼ 0:1 m and Dx ¼ Dy ¼ 0:05 m, were used in the initial computations. Steady state solutions were reached after 20,000 steps. Comparison between numerical results based on different lattice sizes
3532
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
are depicted in Fig. 2, showing that there is little difference. Further quantitative comparison indicates that there is small difference of the results between lattice Dx ¼ 0:2 m and Dx ¼ 0:1 m, but the results are almost the same for both Dx ¼ 0:1 m and Dx ¼ 0:05 m. This suggests that the results based on lattice with Dx 6 0:1 m provide lattice-independent solutions; hence the solution obtained with Dx ¼ 0:05 m is presented here. Fig. 3 shows the profile of the water surface along the channel, in which there is a surface drop because the flow is subcritical. Comparison of the result with the analytical solution is plotted in Fig. 4, showing an
Fig. 2. Subcritical flow: effect of the lattice size on solutions.
Fig. 3. Subcritical flow: profile of the water surface.
Fig. 4. Subcritical flow: comparison of the water surface.
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3533
excellent agreement. The quantitative comparison with the analytical solution indicates that the maximum relative error for the water depth is smaller than 0:25%. In order to test the conservative property of the model, the numerical solution of the discharge (mass) is shown in Fig. 5. Also the quantitative comparison with the theoretical discharge shows that the maximum relative error is smaller than 0:15%. This suggests that the model is accurate and conservative. 3.2. Tidal wave flow Tidal waves often occur in coastal engineering. Here we consider the test problem that Bermudez and V azquez [13] used for verification of an upwind discretisation of the bed slope source terms. This is a one dimensional problem with bed topography defined by (see Fig. 6)
40x 4x 1 H ðxÞ ¼ 50:5 10 sin p ; ð22Þ L L 2
Fig. 5. Subcritical flow: comparison of the discharge.
Fig. 6. Tidal flow: comparison of the water surface at t ¼ 9117:5 s.
3534
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
where L ¼ 14,000 m is the channel length and H ðxÞ the partial depth between a fixed reference level and the bed surface, hence zb ðxÞ ¼ H ð0Þ H ðxÞ. The initial and boundary conditions are
and
hðx; 0Þ ¼ H ðxÞ;
ð23Þ
uðx; 0Þ ¼ 0
ð24Þ
hð0; tÞ ¼ H ð0Þ þ 4 4 sin p
4t 1 þ 86;400 2
;
uðL; tÞ ¼ 0:
ð25Þ ð26Þ
Under these conditions, the tidal wave is relatively short and an asymptotic analytical solution is given by Bermudez and V azquez [13] as
hðx; tÞ ¼ H ðxÞ þ 4 4 sin p
4t 1 þ ; 86;400 2
ðx 14;000Þp 4t 1 cos p uðx; tÞ ¼ þ : 5400hðx; tÞ 86;400 2
ð27Þ ð28Þ
To achieve a lattice-independent solution, three lattices of 600, 800 and 1000 were used. e ¼ 200 m/s and s ¼ 0:6. Eqs. (23)–(26) were used as the initial and boundary conditions. Fig. 7 shows the comparison of the results based on the lattices. Difference in the results with 600 and 800 lattices clearly exists, but the results based on 800 and 1000 lattices are almost the same. Hence the results with 800 lattices are further described here. A comparison of the numerical results with the asymptotic analytical solutions at t ¼ 9117:5 s is shown in Figs. 6 and 7, respectively, showing good agreement. The quantitative comparison with the analytical solution is carried out. For the water depths, the maximum relative error is smaller than 1%. For the velocity, the maximum relative error is smaller than 5% when x 6 11,445 m; when x > 11,445 m, the
Fig. 7. Tidal flow: comparison of the velocity uðx; tÞ at t ¼ 9117:5 s.
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3535
relative error may be larger than 5% and it becomes big with x, which is due to the fact that the zero boundary condition was used for velocity at the downstream end, leading to u ! 0 while x ! 14,000 m. Consequently, the velocity becomes very small such that the relative error is no longer an appropriate measurement for accuracy; instead the absolute error is usually used. For this case, the further calculation shows that the absolute error is smaller than 0:006 m/s for all the velocities with the relative error larger than 5%. Since the analytical solution is obtained by use of the asymptotic analysis, such agreement may be considered to be excellent. This confirms that the model is accurate for unsteady shallow water flow problems. The present method can provide solution of the same accuracy as that reported by Bermudez and Vazquez who used a complex upwind discretisation for the bed slope source terms [13]. 3.3. Flow around a cylinder Flow around a cylinder is a classical problem in fluvial hydraulics. In order to test the capability of the present method, a cylinder with radius of 0.11 m situated in the centre of a channel is used, which is the same as Test 1 used by Yulistiyanto et al. in the numerical and experimental investigations [14]. The channel is 4 m long and 2 m wide (see Fig. 9). The discharge is Q ¼ 0:248 m3 /s; outflow depth is ho ¼ 0:185 m; the bed slope is ozb =ox ¼ 6:25 104 in the flow direction; and the Manning’s coefficient is nb ¼ 0:012. In the computation, 600 300 square lattices were used. Dx ¼ Dy ¼ 0:00667 m, Dt ¼ 0:00145 s and s ¼ 1:982. No-slip boundary condition is used along the wall of the cylinder; slip boundary condition is used for the side walls; the depth at the outflow is set to h ¼ ho ; zero gradient for the depth is specified at the inflow boundary; the inflow velocity u is decided based on the discharge; and v ¼ 0 at the inflow boundary. After 40,000 steps, a steady state solution was reached. Comparison of the depths along the centerline of the channel between the computation and the experimental data is shown in Fig. 8. The calculation indicates that the maximum relative error for the depth is 9.2%; hence agreement is reasonably good. The contours of the water depths are plotted in Fig. 9, showing a similar pattern to Fig. 10 obtained in the experiment. Fig. 11 shows the flow characteristics in the vicinity of the cylinder, which is again similar to Fig. 12 given by Yulistiyanto et al. [14] with a different numerical method. In fact, there is a relation between the wake length Lw and the cylinder diameter D, i.e. Lw 1:3D (see Fig. 12 for definitions). This is also in good agreement with that obtained by Yulistiyanto et al. [14]. The vectors of the whole flow field is depicted in Fig. 13. However, the difference exists in the region close to the cylinder, which may be due to the use of both no-slip boundary conditions and uniform square lattices in the computations. The streamlines are plotted in Fig. 14. 3.4. Flow over a sudden-expansion channel In this section, flow over a channel with a symmetric sudden-expansion is considered to demonstrate that the method is able to simulate a recirculation in shallow water flows. The entrance channel is 1 m wide and
Fig. 8. Flow around cylinder: comparison of the depth along channel centerline.
3536
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
Fig. 9. Flow around cylinder: contours of the water depth (computation).
Fig. 10. Flow around cylinder: contours of the water depth (experiment) from Ref. [14].
Fig. 11. Flow around cylinder: velocity vectors (computation).
Fig. 12. Flow around cylinder: velocity vectors from Ref. [14].
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3537
Fig. 13. Flow around cylinder: velocity vectors.
Fig. 14. Flow around cylinder: streamlines.
Fig. 15. Sudden-expansion channel: velocity vectors.
2 m long, and the expanded channel is 3 m wide and 4 m long (see Fig. 15). There is no bed slope in the channel, and the bed friction is neglected. In the computation, 120 60 square lattices were used. Dx ¼ Dy ¼ 0:05 m, Dt ¼ 0:025 s and s ¼ 1. No-slip boundary condition is used along the walls. The depth at the outflow is set to h ¼ 0:16 m; zero gradient for the depth is specified at the inflow boundary; the discharge is Q ¼ 0:032 m3 /s; and v ¼ 0 at the inflow boundary. After 10,000 steps, a steady state solution was reached. The velocity vectors are plotted in Fig. 15, showing well-developed circulating flows on both sides. This confirms that the present method is capable of simulating recirculation occurring in shallow water flows.
3538
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
4. Conclusions An accurate LABSWE is presented in the paper. The model is simple and conservative. It is able to solve steady and unsteady flow problems. The source terms can naturally be added to the lattice Boltzmann equation as a force term without any special treatment. This makes the model have potential capability in solving practical flow problems.
Acknowledgements We would like to thank Professor W.H. Graf (Laboratoire de Recherches Hydrauliques, Switzerland) for his providing the experimental data of the flow around a cylinder.
Appendix A. 7-Speed model A lattice Boltzmann model for the shallow water equations using the 7-speed hexagonal lattice shown in Fig. A.1 can similarly be developed as 1 Dt fa ðx þ ea Dt; t þ DtÞ fa ðx; tÞ ¼ ðfa faeq Þ þ 2 eai Fi ðx; tÞ s 3e
ðA:1Þ
with the velocity vector of particles defined by ea ¼
0; ~ıe cosðpða1Þ Þ þ~ |e sinðpða1Þ Þ; 3 3
a ¼ 0; a¼16
ðA:2Þ
and the local equilibrium function faeq expressed as
faeq
8 gh2 h > > < h 2 þ 2 ui uj dij ; e e ¼ 2 > gh h 2h h > : þ 2 eai ui þ 4 eai eaj ui uj 2 ui uj dij ; 2 3e 3e 2e 6e
a ¼ 0; ðA:3Þ a ¼ 1 6:
Fig. A.1. 7-Speed hexagonal lattice.
J.G. Zhou / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3527–3539
3539
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]
R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: theory and applications, Phys. Rep. 222 (1992) 45–197. S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Ann. Rev. Fluid Mech. 30 (1998) 329–364. J.G. Zhou, Velocity–depth coupling in shallow water flows, J. Hydr. Engrg., ASCE 121 (10) (1995) 717–724. P.K. Stansby, J.G. Zhou, Shallow-water flow solver with non-hydrostatic pressure: 2D vertical plane problems, Int. J. Numer. Meth. Fluids 28 (1998) 541–563. E.F. Toro, Riemann problems and the WAF method for solving two-dimensional shallow water equations, Phil. Trans. Roy. Soc. Lond. A 338 (1992) 43–68. R.J. LeVeque, Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wavepropagation algorithm, J. Comput. Phys. 146 (1998) 346–365. P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component system, Phys. Rev. 94 (1954) 511–525. P.A. Skordos, Initial and boundary conditions for the lattice Boltzmann method, Phys. Rev. E 48 (1993) 4823–4842. D.H. Rothman, S. Zaleski, Lattice-Gas Cellular Automata, Cambridge University Press, London, 1997. Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Phys. Fluids 9 (1997) 1591– 1598. N. Goutal, F. Maurel (Eds.), Proceedings of the Second Workshop on Dam-break Wave Simulation, HE-43/97/016/B, Departement Laboratoire National d’Hydraulique, Groupe Hydraulique Fluviale, Electricite de France, France, 1997. M.E. V azquez-Cend on, Improved treatment of source terms in upwind schemes for shallow water equations in channels with irregular geometry, J. Comput. Phys. 148 (1999) 497–526. A. Bermudez, M.E. Vazquez, Upwind methods for hyperbolic conservation laws with source terms, Comput. Fluids 23 (1994) 1049–1071. B. Yulistiyanto, Y. Zech, W.H. Graf, Flow around a cylinder: shallow-water modeling with diffusion–dispersion, J. Hydr. Engrg., ASCE 124 (4) (1998) 419–429.