- Email: [email protected]
Enhancement of the LABSWE for shallow water flows
Journal of Computational Physics 230 (2011) 394–401
Contents lists available at ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Enhancement of the LABSWE for shallow water flows Jian Guo Zhou School of Engineering, University of Liverpool, Liverpool L69 3GQ, United Kingdom
a r t i c l e
i n f o
Article history: Received 23 June 2010 Received in revised form 21 September 2010 Accepted 23 September 2010 Available online 25 October 2010 Keywords: Bed slope Force term Lattice Boltzmann method Rectangular lattice Shallow water equations Numerical method Mathematical model
a b s t r a c t In the lattice Boltzmann method for the shallow water equations (LABSWE), the force term involves the first order derivative related to a bed slope, which has to be determined using the centred scheme for an accurate solution. However, such calculation contradicts the spirit of only simple arithmetic calculations required in the lattice Boltzmann method. This paper shows how to remove this drawback from the LABSWE by incorporating the bed level into the lattice Boltzmann equation in a consistent manner with the lattice Boltzmann approach. Three numerical tests have been presented to verify the proposed method. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction Shallow water equations describe a wide range of fluid flows, e.g. tidal flows, waves, open channel flows, dam breaks, and atmospheric flows. There are many numerical methods available for solution of the equations, which turns out to be a very successful tool in studying diverse flow problems encountered in engineering [1–7]. As natural shallow flows always occur over complex topography, an accurate treatment of the bed slope becomes the key to success or failure in a numerical method. For instance, in a Godunov-type method, a special treatment for bed slope has to be applied for accurate solutions [8–10]; in the lattice Boltzmann method for shallow water equations (LABSWE), the use of the centred scheme for calculation of a bed slope is necessary to produce accurate solutions [7]. As an alternative computational method, the LABSWE has been applied to several complex flow problems, demonstrating its potential, capability and accuracy in simulating shallow water flows [11–13]. However, in the LABSWE, the centred scheme involves the determination of the first order derivative associated with a bed slope, which is inconsistent with the spirit of the lattice Boltzmann hydrodynamics. In this paper, a novel incorporation of the bed level into the lattice Boltzmann equation is proposed to eliminate the calculation of the derivative, hence enhancing the LABSWE for predicting shallow water flows. Three numerical tests are carried out and compared with either analytical solutions or other numerical results to validate the accuracy and capability of the method. 2. Improved LABSWE In the LABSWE [7], the lattice Boltzmann equation including a force term on a square lattice with nine particle velocities (D2Q9) shown in Fig. 1 is
E-mail address: [email protected] 0021-9991/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2010.09.027
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
4
3
5 6
395
2 1
7
8
Fig. 1. Nine-velocity square lattice (D2Q9).
fa ðx þ ea Dt; t þ DtÞ fa ðx; tÞ ¼
1
s
Dt faeq fa þ 2 eaj F j ; 6e
ð1Þ
where fa is the particle distribution function; x is the space vector defined by Cartesian coordinates, i.e., x = (x, y) in 2D space; t is the time; Dt is the time step; ea is the particle velocity vector; eaj is the component of ea in j direction; e = Dx/Dt, Dx is the lattice size; s is the single relaxation time [14]; and faeq is the local equilibrium distribution function defined as
faeq
8 > i ui < h 1 5gh 2u a ¼ 0; ; 3e2 6e2 ¼ e e u u > : ka h gh2 þ eai u2 i þ ai aj4 i j ui u2i ; a – 0; 6e 3e 2e 6e
ð2Þ
in which ka = 1 when a = 1, 3, 5, 7 and ka = 1/4 when a = 2, 4, 6, 8; and Fi is the force term given by,
F i ¼ gh
@zb swi sbi þ þ Xhuy dix Xhux diy ; @xi q q
ð3Þ
where, h and ui are the respective water depth and the depth-averaged velocity; g = 9.81 m/s2 is the gravitational acceleration; zb is the bed elevation above datum; swi is the wind shear stress; sbi is the bed shear stress; q is the water density; X is the Coriolis parameter for the effect of the earth’s rotation; and dij is the Kronecker delta function,
dij ¼
0;
i – j;
ð4Þ
1; i ¼ j:
The physical variables of water depth and velocity can be calculated as
h¼
X
fa ;
ð5Þ
a
and
ui ¼
1X eai fa : h a
ð6Þ
It has been shown from the Chapman-Enskog analysis that the LABSWE can produce an accurate solution to the shallow water equations,
@h @ðhuj Þ ¼ 0; þ @t @xj
ð7Þ 2
@ðhui Þ @ðhui uj Þ g @h @ 2 ðhui Þ þ ¼ þm þ Fi; @t @xj 2 @xi @x2j
ð8Þ
where i and j are indices; xi is the Cartesian coordinate; and m is the kinematic viscosity. b As seen from Eq. (3), Fi contains a first derivative related to the bed slope, gh @z , which makes the model inconsistent with @xi the lattice Boltzmann hydrodynamics. To remedy this, we propose a lattice Boltzmann equation,
fa ðx þ ea Dt; t þ DtÞ fa ðx; tÞ ¼
1
s
gh Dt faeq fa 2 ½zb ðx þ ea DtÞ zb ðxÞ þ 2 eaj F j ; 6e 6e
ð9Þ
396
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
¼ 0:5½hðx þ ea Dt; t þ DtÞ þ hðx; tÞ, to simulate the following shallow water equations in which h
@h @ðhuj Þ þ ¼ 0; @t @xj
ð10Þ 2
@ðhui Þ @ðhui uj Þ g @h @zb @ 2 ðhui Þ þ ¼ gh þm þ Fi; @t @xj 2 @xi @xi @x2j
ð11Þ
with Fi redefined as
Fi ¼
swi sbi þ Xhuy dix Xhux diy : q q
ð12Þ
By using the Chapman-Enskog analysis, it can be proved that the shallow water Eqs. (10) and (11) can be recovered from the lattice Boltzmann Eq. (9). Assuming that Dt is small and
Dt ¼ e;
ð13Þ
we have Eq. (9) expressed as
1 gh e fa ðx þ ea e; t þ eÞ fa ðx; tÞ ¼ ðfa faeq Þ 2 ½zb ðx þ ea eÞ zb ðxÞ þ 2 eaj F j : 6e s 6e
ð14Þ
Taking a Taylor expansion to the first term on the left-hand side of the above equation in time and space around point (x, t) leads to
2 @ @ e2 @ @ fa þ þ eaj þ eaj fa þ Oðe3 Þ: @t @xj @xj 2 @t
e
ð15Þ ð0Þ
According to Chapman-Enskog procedure, fa can be expanded around fa ,
fa ¼ fað0Þ þ efað1Þ þ e2 fað2Þ þ Oðe2 Þ:
ð16Þ
The second term on the right hand side of Eq. (14) can also be expressed via the Taylor expansion,
g e @h @h þ e h þ a j 6e2 @xj 2 @t
@z e2 @ 2 zb eeaj b þ eai eaj @xj 2 @xi @xj
! þ Oðe3 Þ:
ð17Þ
The centred scheme [7] is used for force term Fj,
1 1 F j ¼ F j x þ ea e; t þ e ; 2 2
ð18Þ
which can again be written, via a Taylor expansion, as
1 1 e @F j @F j þ Oðe2 Þ: F j x þ ea e; t þ e ¼ F j þ þ eai 2 2 2 @t @xi
ð19Þ
After inserting Eqs. (15), (16), (17) and (19) into Eq. (14), the equation to order e0 is
fað0Þ ¼ faeq ;
ð20Þ
to order e
ð1Þ gheaj @zb eaj F j @ @ fa fað0Þ ¼ þ 2 ; þ eaj @t @xj s 6e2 @xj 6e
ð21Þ
and to order e2 as
2 geaj @h @ @ 1 @ @ 1 @h @zb gheai eaj @ 2 zb eaj @F j @F j fað1Þ þ : þ eaj þ eaj þ e fað0Þ ¼ fað2Þ þ þ e a i a i @t @xj 2 @t @xj @xi @xj s 12e2 @t 12e2 @xi @xj 12e2 @t @xi ð22Þ
Substitution of Eq. (21) into Eq. (22) gives
1 @ @ 1 f ð1Þ ¼ fað2Þ : þ eaj 1 2s @t @xj a s Taking
P
ð23Þ
[(21) + e (23)] about a provides
@ X ð0Þ @ X f þ eaj fað0Þ ¼ 0: @t a a @xj a
ð24Þ
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
397
Evaluation of the terms in the above equation using Eq. (2) results in the second-order accurate continuity Eq. (10). P Taking eai [(21) + e (23)] about a yields
@ X @ X 1 @ X @zb eai fað0Þ þ eai eaj fað0Þ þ e 1 eai eaj fað1Þ ¼ gh þ Fi: @t a @xj a 2s @xj a @xi
ð25Þ
After the terms are evaluated with Eq. (2) and some algebra, the above equation becomes the momentum equation (11), which is second-order accurate. can be eliminated by using the method by He et al. [15]; (b) It must be pointed out that (a) the implicitness related to h alternatively, the following semi-implicit form,
¼ 0:5½hðx þ ea Dt; tÞ þ hðx; tÞ; h
ð26Þ
can be used, which is simple and demonstrated to produce accurate solutions, and hence it is preferred in practice; and (c) for the rectangular lattice Boltzmann method [16], Eq. (9) should be replaced by
fa ðx þ ea Dt; t þ DtÞ fa ðx; tÞ ¼
1
s
ðfaeq fa Þ X a þ Z a Dt;
ð27Þ
in which
Xa ¼
8 gh > a ¼ 1; 5; > 2 ½zb ðx þ ea DtÞ zb ðxÞ; > < 2ex
gh ½z ðx þ ea DtÞ zb ðxÞ; a ¼ 3; 7; 2e2y b > > > : 0; otherwise;
ð28Þ
and
Za ¼
8 0; > > > Fx > > < 6eax ;
a ¼ 0; a ¼ 1 & 5; ð29Þ
F
y ; a ¼ 3 & 7; > 6eay > > > > : F i ; otherwise; 6e
ai
where ex = Dx/Dt and ey = Dy/Dt with lattice sizes of Dx and Dy in x and y directions, respectively. It is worthy of mention that both the present method and the centred scheme [7] can provide an accurate solution up to computer precision, and hence the results are effectively identical as demonstrated in the following numerical validation. However, the advantage of the proposed method over the centred scheme lies in the fact that it completely retains the simple arithmetic operations in the LABSWE without calculation of a derivative, leading to easier implementation in the algorithm.
3. Validation In order to verify the described model, three numerical tests are presented. The SI units are used for the physical variables in the numerical simulations. 3.1. Stationary case First of all, a 1D stationary case is investigated. As there is a non-zero force associated with the bed slope, many numerical methods fail to reproduce the exact solution and it becomes a benchmark test to validate a numerical method for shallow water flows; hence it is chosen as the first test. The bed topography is the same as that used by LeVeque [9] and is given by
( h zb ðxÞ ¼
1 4
0;
i cos pðx0:5Þ þ 1 ; jx 0:5j < 0:1; 0:1
ð30Þ
otherwise:
The channel has the length of 1. In the numerical computation, 400 lattices are used: Dx = 0.0025, Dt = 0.0005 and s = 0.8. The initial conditions are ui = 0 and h + zb = 1, for which the exact solution is ui 0 and h + zb 1. This is a severe benchmark test because even a little inappropriate formulation in a numerical method can generate an artificial flow [8]. Hence the case is an appropriate test for the proposed scheme. During the simulation, it is found that the accurate results up to computer precision for both water depth and velocity have been retained for any running period since the start of the computation, indicating that the proposed scheme is accurate. Comparisons of water level and velocity with the exact solutions are plotted in Figs. 2 and 3, showing perfect agreements.
398
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
Fig. 2. Comparison of water level.
Fig. 3. Comparison of velocity.
3.2. Tidal flow Then, a tidal flow over an irregular bed is predicted, which is a common flow problem in coastal engineering. The bed is defined with data listed in Table 1. This is also a 1D problem with the initial and boundary conditions given by
hðx; 0Þ ¼ 16 zb ðxÞ;
ð31Þ
ux ðx; 0Þ ¼ 0;
ð32Þ
and
hð0; tÞ ¼ 20 4 sin
p
4t 1 þ 86; 400 2
ð33Þ
;
ux ð1500; tÞ ¼ 0:
ð34Þ
In the simulation, 200 lattices are used with Dx = 3.75, Dt = 0.15 and s = 1.5. Two numerical results at t = 10,800 and t = 32,400 corresponding to the half-risen tidal flow with maximum positive velocities and to the half-ebb tidal flow with maximum negative velocities are compared with the analytical solutions [8] and depicted in Figs. 4 and 5, respectively. The maximum relative errors are less than 0.005% for the water level, less than 0.05% for velocity larger than 0.002, and less than 0.3% for smaller velocity, revealing excellent agreements. It is checked that the LABSWE with the centred scheme [7]
Table 1 Bed elevation zb for irregular bed. x (m) zb (m)
0 0
50 0
100 2.5
150 5
250 5
300 3
350 5
400 5
425 7.5
435 8
x (m) zb (m)
450 9
475 9
500 9.1
505 9
530 9
550 6
565 5.5
575 5.5
600 5
650 4
x (m) zb (m)
700 3
750 3
800 2.3
820 2
900 1.2
950 0.4
1000 0
1500 0
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
399
Fig. 4. Comparison of velocity at t = 10,800 s.
Fig. 5. Comparison of velocity at t = 32,400 s.
can also provide the similar accurate results, which is a natural proof that both the present method and the centred scheme can generate accurate solutions. 3.3. Wind-driven circulation Finally, we consider a wind-driven circulation in a lake, which may generate a complex flow phenomenon depending on the bed topography of a lake. In this test, a uniform wind shear stress is applied to the shallow water in a circular basin with the bed topography defined by the still water depth H,
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 x2 þ y2 A 1 @1 1 ; þ Hðx; yÞ ¼ 386:4 1:3 2 2
ð35Þ
from which, the bed level can be determined as zb(x, y) = H(0, 0) H(x, y). The same dish-shaped basin is also used by Rogers et al. [17] to test a Godunov-type method. Initially, the water in the basin is still and then a uniform wind shear stress from southwest to northeast is applied, i.e.
p p ðswx ; swy Þ ¼ 0:004q cos ; sin : 4 4
ð36Þ
As a result, the steady flow consists of two relatively strong counter-rotating gyres with flow in the deeper water against the direction of the wind. This exhibits complex flow phenomenon. In the numerical computation, 200 200 lattices are used. Dx = 2, Dt = 0.2 and s = 1.325. After the steady solution is obtained, the flow field is shown in Fig. 6 and the normalised resultant velocities at cross section A–A are compared with the analytical solution [18] in Fig. 7, exhibiting similar agreement to that by Rogers et al. using a high-resolution Godunov-type method [17]. Although there is discrepancy between the numerical prediction and the analytical solution, such agreement is reasonable due to the fact that the assumptions of both the rigid-lid approximation for the water surface and a parabolic distribution for the eddy viscosity were used in deriving the analytical solution. Furthermore, the original numerical results [7] from the LABSWE with the centred scheme for the bed slope are also plotted in the figure, showing very good agreement. In fact, both numerical results are almost identical except for the small region close to the boundary where the maximum relative error is less than 0.1%, confirming again the accuracy and correctness of the proposed scheme.
400
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
Fig. 6. Flow field for wind-driven flow.
Fig. 7. Comparison of the resultant velocities along Cross-section A–A (see Fig. 6) and U0 = 0.89, s = ux + uy.
4. Conclusions The paper presents a novel treatment for the first order derivative associated with the bed slope in the LABSWE. The scheme preserves the simple arithmetic calculations of the lattice Boltzmann method. Numerical tests have shown that the method can provide accurate solutions, making the LABSWE an ideal model for simulating shallow water flows. References [1] A. Alcrudo, P. Garcia-Navarro, A high resolution Godunov-type scheme in finite volumes for the 2D shallow water equations, International Journal for Numerical Methods in Fluids 16 (1993) 489–505. [2] V. Casulli, Semi-implicit finite difference methods for the two-dimensional shallow water equations, Journal of Computational Physics 86 (1990) 56– 74. [3] C.B. Vreugdenhil, Numerical Methods for Shallow-Water Flow, Kluwer Academic Publishers, Dordrecht, 1994. [4] A.G.L. Borthwick, G.A. Akponasa, Reservoir flow prediction by contravariant shallow water equations, Journal of Hydraulic Engineering, ASCE 123 (5) (1997) 432–439. [5] B. Yulistiyanto, Y. Zech, W.H. Graf, Flow around a cylinder: shallow-water modeling with diffusion-dispersion, Journal of Hydraulic Engineering, ASCE 124 (4) (1998) 419–429. [6] K. Hu, C.G. Mingham, D.M. Causon, Numerical simulation of wave overtopping of coastal structures using the non-linear shallow water equations, Coastal Engineering 41 (4) (2000) 433–465. [7] J.G. Zhou, Lattice Boltzmann Methods for Shallow Water Flows, Springer-Verlag, Berlin, 2004. [8] A. Bermudez, M.E. Vázquez, Upwind methods for hyperbolic conservation laws with source terms, Computers and Fluids 23 (1994) 1049–1071. [9] R.J. LeVeque, Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm, Journal of Computational Physics 146 (1998) 346–365.
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
401
[10] J.G. Zhou, D.M. Causon, C.G. Mingham, D.M. Ingram, The surface gradient method for the treatment of source terms in the shallow-water equations, Journal of Computational Physics 168 (2001) 1–25. [11] H. Liu, J.G. Zhou, R. Burrows, Lattice boltzmann model for shallow water flows in curved and meandering channels, International Journal of Computational Fluid Dynamics 23 (2009) 209–220. [12] H. Liu, J.G. Zhou, R. Burrows, Multi-block lattice boltzmann simulations of subcritical flow in open channel junctions, Computers and Fluids 38 (2009) 1108–1117. [13] H. Liu, J.G. Zhou, R. Burrows, Numerical modelling of turbulent compound channel flow using the lattice boltzmann method, International Journal for Numerical Methods in Fluids 59 (2009) 753–765. [14] 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, Physical Review 94 (3) (1954) 511–525. [15] X. He, S. Chen, G.D. Doolen, A novel thermal model for the lattice boltzmann method in incompressible limit, Journal of Computational Physics 146 (1998) 282–300. [16] J.G. Zhou, Rectangular lattice Boltzmann method, Physical Review E 81 (2010) 026705. [17] B. Rogers, M. Fujihara, A.G.L. Borthwick, Adaptive Q-tree Godunov-type scheme for shallow water equations, International Journal for Numerical Methods in Fluids 35 (2001) 247–280. [18] C. Kranenburg, Wind-driven chaotic advection in a shallow model lake, Journal of Hydraulic Research 30 (1) (1992) 29–46.
Contents lists available at ScienceDirect
Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp
Enhancement of the LABSWE for shallow water flows Jian Guo Zhou School of Engineering, University of Liverpool, Liverpool L69 3GQ, United Kingdom
a r t i c l e
i n f o
Article history: Received 23 June 2010 Received in revised form 21 September 2010 Accepted 23 September 2010 Available online 25 October 2010 Keywords: Bed slope Force term Lattice Boltzmann method Rectangular lattice Shallow water equations Numerical method Mathematical model
a b s t r a c t In the lattice Boltzmann method for the shallow water equations (LABSWE), the force term involves the first order derivative related to a bed slope, which has to be determined using the centred scheme for an accurate solution. However, such calculation contradicts the spirit of only simple arithmetic calculations required in the lattice Boltzmann method. This paper shows how to remove this drawback from the LABSWE by incorporating the bed level into the lattice Boltzmann equation in a consistent manner with the lattice Boltzmann approach. Three numerical tests have been presented to verify the proposed method. Ó 2010 Elsevier Inc. All rights reserved.
1. Introduction Shallow water equations describe a wide range of fluid flows, e.g. tidal flows, waves, open channel flows, dam breaks, and atmospheric flows. There are many numerical methods available for solution of the equations, which turns out to be a very successful tool in studying diverse flow problems encountered in engineering [1–7]. As natural shallow flows always occur over complex topography, an accurate treatment of the bed slope becomes the key to success or failure in a numerical method. For instance, in a Godunov-type method, a special treatment for bed slope has to be applied for accurate solutions [8–10]; in the lattice Boltzmann method for shallow water equations (LABSWE), the use of the centred scheme for calculation of a bed slope is necessary to produce accurate solutions [7]. As an alternative computational method, the LABSWE has been applied to several complex flow problems, demonstrating its potential, capability and accuracy in simulating shallow water flows [11–13]. However, in the LABSWE, the centred scheme involves the determination of the first order derivative associated with a bed slope, which is inconsistent with the spirit of the lattice Boltzmann hydrodynamics. In this paper, a novel incorporation of the bed level into the lattice Boltzmann equation is proposed to eliminate the calculation of the derivative, hence enhancing the LABSWE for predicting shallow water flows. Three numerical tests are carried out and compared with either analytical solutions or other numerical results to validate the accuracy and capability of the method. 2. Improved LABSWE In the LABSWE [7], the lattice Boltzmann equation including a force term on a square lattice with nine particle velocities (D2Q9) shown in Fig. 1 is
E-mail address: [email protected] 0021-9991/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2010.09.027
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
4
3
5 6
395
2 1
7
8
Fig. 1. Nine-velocity square lattice (D2Q9).
fa ðx þ ea Dt; t þ DtÞ fa ðx; tÞ ¼
1
s
Dt faeq fa þ 2 eaj F j ; 6e
ð1Þ
where fa is the particle distribution function; x is the space vector defined by Cartesian coordinates, i.e., x = (x, y) in 2D space; t is the time; Dt is the time step; ea is the particle velocity vector; eaj is the component of ea in j direction; e = Dx/Dt, Dx is the lattice size; s is the single relaxation time [14]; and faeq is the local equilibrium distribution function defined as
faeq
8 > i ui < h 1 5gh 2u a ¼ 0; ; 3e2 6e2 ¼ e e u u > : ka h gh2 þ eai u2 i þ ai aj4 i j ui u2i ; a – 0; 6e 3e 2e 6e
ð2Þ
in which ka = 1 when a = 1, 3, 5, 7 and ka = 1/4 when a = 2, 4, 6, 8; and Fi is the force term given by,
F i ¼ gh
@zb swi sbi þ þ Xhuy dix Xhux diy ; @xi q q
ð3Þ
where, h and ui are the respective water depth and the depth-averaged velocity; g = 9.81 m/s2 is the gravitational acceleration; zb is the bed elevation above datum; swi is the wind shear stress; sbi is the bed shear stress; q is the water density; X is the Coriolis parameter for the effect of the earth’s rotation; and dij is the Kronecker delta function,
dij ¼
0;
i – j;
ð4Þ
1; i ¼ j:
The physical variables of water depth and velocity can be calculated as
h¼
X
fa ;
ð5Þ
a
and
ui ¼
1X eai fa : h a
ð6Þ
It has been shown from the Chapman-Enskog analysis that the LABSWE can produce an accurate solution to the shallow water equations,
@h @ðhuj Þ ¼ 0; þ @t @xj
ð7Þ 2
@ðhui Þ @ðhui uj Þ g @h @ 2 ðhui Þ þ ¼ þm þ Fi; @t @xj 2 @xi @x2j
ð8Þ
where i and j are indices; xi is the Cartesian coordinate; and m is the kinematic viscosity. b As seen from Eq. (3), Fi contains a first derivative related to the bed slope, gh @z , which makes the model inconsistent with @xi the lattice Boltzmann hydrodynamics. To remedy this, we propose a lattice Boltzmann equation,
fa ðx þ ea Dt; t þ DtÞ fa ðx; tÞ ¼
1
s
gh Dt faeq fa 2 ½zb ðx þ ea DtÞ zb ðxÞ þ 2 eaj F j ; 6e 6e
ð9Þ
396
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
¼ 0:5½hðx þ ea Dt; t þ DtÞ þ hðx; tÞ, to simulate the following shallow water equations in which h
@h @ðhuj Þ þ ¼ 0; @t @xj
ð10Þ 2
@ðhui Þ @ðhui uj Þ g @h @zb @ 2 ðhui Þ þ ¼ gh þm þ Fi; @t @xj 2 @xi @xi @x2j
ð11Þ
with Fi redefined as
Fi ¼
swi sbi þ Xhuy dix Xhux diy : q q
ð12Þ
By using the Chapman-Enskog analysis, it can be proved that the shallow water Eqs. (10) and (11) can be recovered from the lattice Boltzmann Eq. (9). Assuming that Dt is small and
Dt ¼ e;
ð13Þ
we have Eq. (9) expressed as
1 gh e fa ðx þ ea e; t þ eÞ fa ðx; tÞ ¼ ðfa faeq Þ 2 ½zb ðx þ ea eÞ zb ðxÞ þ 2 eaj F j : 6e s 6e
ð14Þ
Taking a Taylor expansion to the first term on the left-hand side of the above equation in time and space around point (x, t) leads to
2 @ @ e2 @ @ fa þ þ eaj þ eaj fa þ Oðe3 Þ: @t @xj @xj 2 @t
e
ð15Þ ð0Þ
According to Chapman-Enskog procedure, fa can be expanded around fa ,
fa ¼ fað0Þ þ efað1Þ þ e2 fað2Þ þ Oðe2 Þ:
ð16Þ
The second term on the right hand side of Eq. (14) can also be expressed via the Taylor expansion,
g e @h @h þ e h þ a j 6e2 @xj 2 @t
@z e2 @ 2 zb eeaj b þ eai eaj @xj 2 @xi @xj
! þ Oðe3 Þ:
ð17Þ
The centred scheme [7] is used for force term Fj,
1 1 F j ¼ F j x þ ea e; t þ e ; 2 2
ð18Þ
which can again be written, via a Taylor expansion, as
1 1 e @F j @F j þ Oðe2 Þ: F j x þ ea e; t þ e ¼ F j þ þ eai 2 2 2 @t @xi
ð19Þ
After inserting Eqs. (15), (16), (17) and (19) into Eq. (14), the equation to order e0 is
fað0Þ ¼ faeq ;
ð20Þ
to order e
ð1Þ gheaj @zb eaj F j @ @ fa fað0Þ ¼ þ 2 ; þ eaj @t @xj s 6e2 @xj 6e
ð21Þ
and to order e2 as
2 geaj @h @ @ 1 @ @ 1 @h @zb gheai eaj @ 2 zb eaj @F j @F j fað1Þ þ : þ eaj þ eaj þ e fað0Þ ¼ fað2Þ þ þ e a i a i @t @xj 2 @t @xj @xi @xj s 12e2 @t 12e2 @xi @xj 12e2 @t @xi ð22Þ
Substitution of Eq. (21) into Eq. (22) gives
1 @ @ 1 f ð1Þ ¼ fað2Þ : þ eaj 1 2s @t @xj a s Taking
P
ð23Þ
[(21) + e (23)] about a provides
@ X ð0Þ @ X f þ eaj fað0Þ ¼ 0: @t a a @xj a
ð24Þ
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
397
Evaluation of the terms in the above equation using Eq. (2) results in the second-order accurate continuity Eq. (10). P Taking eai [(21) + e (23)] about a yields
@ X @ X 1 @ X @zb eai fað0Þ þ eai eaj fað0Þ þ e 1 eai eaj fað1Þ ¼ gh þ Fi: @t a @xj a 2s @xj a @xi
ð25Þ
After the terms are evaluated with Eq. (2) and some algebra, the above equation becomes the momentum equation (11), which is second-order accurate. can be eliminated by using the method by He et al. [15]; (b) It must be pointed out that (a) the implicitness related to h alternatively, the following semi-implicit form,
¼ 0:5½hðx þ ea Dt; tÞ þ hðx; tÞ; h
ð26Þ
can be used, which is simple and demonstrated to produce accurate solutions, and hence it is preferred in practice; and (c) for the rectangular lattice Boltzmann method [16], Eq. (9) should be replaced by
fa ðx þ ea Dt; t þ DtÞ fa ðx; tÞ ¼
1
s
ðfaeq fa Þ X a þ Z a Dt;
ð27Þ
in which
Xa ¼
8 gh > a ¼ 1; 5; > 2 ½zb ðx þ ea DtÞ zb ðxÞ; > < 2ex
gh ½z ðx þ ea DtÞ zb ðxÞ; a ¼ 3; 7; 2e2y b > > > : 0; otherwise;
ð28Þ
and
Za ¼
8 0; > > > Fx > > < 6eax ;
a ¼ 0; a ¼ 1 & 5; ð29Þ
F
y ; a ¼ 3 & 7; > 6eay > > > > : F i ; otherwise; 6e
ai
where ex = Dx/Dt and ey = Dy/Dt with lattice sizes of Dx and Dy in x and y directions, respectively. It is worthy of mention that both the present method and the centred scheme [7] can provide an accurate solution up to computer precision, and hence the results are effectively identical as demonstrated in the following numerical validation. However, the advantage of the proposed method over the centred scheme lies in the fact that it completely retains the simple arithmetic operations in the LABSWE without calculation of a derivative, leading to easier implementation in the algorithm.
3. Validation In order to verify the described model, three numerical tests are presented. The SI units are used for the physical variables in the numerical simulations. 3.1. Stationary case First of all, a 1D stationary case is investigated. As there is a non-zero force associated with the bed slope, many numerical methods fail to reproduce the exact solution and it becomes a benchmark test to validate a numerical method for shallow water flows; hence it is chosen as the first test. The bed topography is the same as that used by LeVeque [9] and is given by
( h zb ðxÞ ¼
1 4
0;
i cos pðx0:5Þ þ 1 ; jx 0:5j < 0:1; 0:1
ð30Þ
otherwise:
The channel has the length of 1. In the numerical computation, 400 lattices are used: Dx = 0.0025, Dt = 0.0005 and s = 0.8. The initial conditions are ui = 0 and h + zb = 1, for which the exact solution is ui 0 and h + zb 1. This is a severe benchmark test because even a little inappropriate formulation in a numerical method can generate an artificial flow [8]. Hence the case is an appropriate test for the proposed scheme. During the simulation, it is found that the accurate results up to computer precision for both water depth and velocity have been retained for any running period since the start of the computation, indicating that the proposed scheme is accurate. Comparisons of water level and velocity with the exact solutions are plotted in Figs. 2 and 3, showing perfect agreements.
398
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
Fig. 2. Comparison of water level.
Fig. 3. Comparison of velocity.
3.2. Tidal flow Then, a tidal flow over an irregular bed is predicted, which is a common flow problem in coastal engineering. The bed is defined with data listed in Table 1. This is also a 1D problem with the initial and boundary conditions given by
hðx; 0Þ ¼ 16 zb ðxÞ;
ð31Þ
ux ðx; 0Þ ¼ 0;
ð32Þ
and
hð0; tÞ ¼ 20 4 sin
p
4t 1 þ 86; 400 2
ð33Þ
;
ux ð1500; tÞ ¼ 0:
ð34Þ
In the simulation, 200 lattices are used with Dx = 3.75, Dt = 0.15 and s = 1.5. Two numerical results at t = 10,800 and t = 32,400 corresponding to the half-risen tidal flow with maximum positive velocities and to the half-ebb tidal flow with maximum negative velocities are compared with the analytical solutions [8] and depicted in Figs. 4 and 5, respectively. The maximum relative errors are less than 0.005% for the water level, less than 0.05% for velocity larger than 0.002, and less than 0.3% for smaller velocity, revealing excellent agreements. It is checked that the LABSWE with the centred scheme [7]
Table 1 Bed elevation zb for irregular bed. x (m) zb (m)
0 0
50 0
100 2.5
150 5
250 5
300 3
350 5
400 5
425 7.5
435 8
x (m) zb (m)
450 9
475 9
500 9.1
505 9
530 9
550 6
565 5.5
575 5.5
600 5
650 4
x (m) zb (m)
700 3
750 3
800 2.3
820 2
900 1.2
950 0.4
1000 0
1500 0
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
399
Fig. 4. Comparison of velocity at t = 10,800 s.
Fig. 5. Comparison of velocity at t = 32,400 s.
can also provide the similar accurate results, which is a natural proof that both the present method and the centred scheme can generate accurate solutions. 3.3. Wind-driven circulation Finally, we consider a wind-driven circulation in a lake, which may generate a complex flow phenomenon depending on the bed topography of a lake. In this test, a uniform wind shear stress is applied to the shallow water in a circular basin with the bed topography defined by the still water depth H,
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 x2 þ y2 A 1 @1 1 ; þ Hðx; yÞ ¼ 386:4 1:3 2 2
ð35Þ
from which, the bed level can be determined as zb(x, y) = H(0, 0) H(x, y). The same dish-shaped basin is also used by Rogers et al. [17] to test a Godunov-type method. Initially, the water in the basin is still and then a uniform wind shear stress from southwest to northeast is applied, i.e.
p p ðswx ; swy Þ ¼ 0:004q cos ; sin : 4 4
ð36Þ
As a result, the steady flow consists of two relatively strong counter-rotating gyres with flow in the deeper water against the direction of the wind. This exhibits complex flow phenomenon. In the numerical computation, 200 200 lattices are used. Dx = 2, Dt = 0.2 and s = 1.325. After the steady solution is obtained, the flow field is shown in Fig. 6 and the normalised resultant velocities at cross section A–A are compared with the analytical solution [18] in Fig. 7, exhibiting similar agreement to that by Rogers et al. using a high-resolution Godunov-type method [17]. Although there is discrepancy between the numerical prediction and the analytical solution, such agreement is reasonable due to the fact that the assumptions of both the rigid-lid approximation for the water surface and a parabolic distribution for the eddy viscosity were used in deriving the analytical solution. Furthermore, the original numerical results [7] from the LABSWE with the centred scheme for the bed slope are also plotted in the figure, showing very good agreement. In fact, both numerical results are almost identical except for the small region close to the boundary where the maximum relative error is less than 0.1%, confirming again the accuracy and correctness of the proposed scheme.
400
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
Fig. 6. Flow field for wind-driven flow.
Fig. 7. Comparison of the resultant velocities along Cross-section A–A (see Fig. 6) and U0 = 0.89, s = ux + uy.
4. Conclusions The paper presents a novel treatment for the first order derivative associated with the bed slope in the LABSWE. The scheme preserves the simple arithmetic calculations of the lattice Boltzmann method. Numerical tests have shown that the method can provide accurate solutions, making the LABSWE an ideal model for simulating shallow water flows. References [1] A. Alcrudo, P. Garcia-Navarro, A high resolution Godunov-type scheme in finite volumes for the 2D shallow water equations, International Journal for Numerical Methods in Fluids 16 (1993) 489–505. [2] V. Casulli, Semi-implicit finite difference methods for the two-dimensional shallow water equations, Journal of Computational Physics 86 (1990) 56– 74. [3] C.B. Vreugdenhil, Numerical Methods for Shallow-Water Flow, Kluwer Academic Publishers, Dordrecht, 1994. [4] A.G.L. Borthwick, G.A. Akponasa, Reservoir flow prediction by contravariant shallow water equations, Journal of Hydraulic Engineering, ASCE 123 (5) (1997) 432–439. [5] B. Yulistiyanto, Y. Zech, W.H. Graf, Flow around a cylinder: shallow-water modeling with diffusion-dispersion, Journal of Hydraulic Engineering, ASCE 124 (4) (1998) 419–429. [6] K. Hu, C.G. Mingham, D.M. Causon, Numerical simulation of wave overtopping of coastal structures using the non-linear shallow water equations, Coastal Engineering 41 (4) (2000) 433–465. [7] J.G. Zhou, Lattice Boltzmann Methods for Shallow Water Flows, Springer-Verlag, Berlin, 2004. [8] A. Bermudez, M.E. Vázquez, Upwind methods for hyperbolic conservation laws with source terms, Computers and Fluids 23 (1994) 1049–1071. [9] R.J. LeVeque, Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm, Journal of Computational Physics 146 (1998) 346–365.
J.G. Zhou / Journal of Computational Physics 230 (2011) 394–401
401
[10] J.G. Zhou, D.M. Causon, C.G. Mingham, D.M. Ingram, The surface gradient method for the treatment of source terms in the shallow-water equations, Journal of Computational Physics 168 (2001) 1–25. [11] H. Liu, J.G. Zhou, R. Burrows, Lattice boltzmann model for shallow water flows in curved and meandering channels, International Journal of Computational Fluid Dynamics 23 (2009) 209–220. [12] H. Liu, J.G. Zhou, R. Burrows, Multi-block lattice boltzmann simulations of subcritical flow in open channel junctions, Computers and Fluids 38 (2009) 1108–1117. [13] H. Liu, J.G. Zhou, R. Burrows, Numerical modelling of turbulent compound channel flow using the lattice boltzmann method, International Journal for Numerical Methods in Fluids 59 (2009) 753–765. [14] 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, Physical Review 94 (3) (1954) 511–525. [15] X. He, S. Chen, G.D. Doolen, A novel thermal model for the lattice boltzmann method in incompressible limit, Journal of Computational Physics 146 (1998) 282–300. [16] J.G. Zhou, Rectangular lattice Boltzmann method, Physical Review E 81 (2010) 026705. [17] B. Rogers, M. Fujihara, A.G.L. Borthwick, Adaptive Q-tree Godunov-type scheme for shallow water equations, International Journal for Numerical Methods in Fluids 35 (2001) 247–280. [18] C. Kranenburg, Wind-driven chaotic advection in a shallow model lake, Journal of Hydraulic Research 30 (1) (1992) 29–46.