A multispeed Discrete Boltzmann Model for ... - Semantic Scholar
Journal of Computational Physics 284 (2015) 117–132
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
A multispeed Discrete Boltzmann Model for transcritical 2D shallow water flows Michele La Rocca a,∗ , Andrea Montessori a , Pietro Prestininzi a , Sauro Succi b a b
Dipartimento di Ingegneria, Universita’ degli Studi Roma TRE, Via Vito Volterra 62, 00146 Rome, Italy Istituto per le Applicazioni del Calcolo, CNR, Via dei Taurini 19, 00189 Rome, Italy
a r t i c l e
i n f o
Article history: Received 23 June 2014 Received in revised form 29 October 2014 Accepted 16 December 2014 Available online 23 December 2014 Keywords: Shallow water equations Multispeed Discrete Boltzmann Model Transcritical flows
a b s t r a c t In this work a Discrete Boltzmann Model for the solution of transcritical 2D shallow water flows is presented and validated. In order to provide the model with transcritical capabilities, a particular multispeed velocity set has been employed for the discretization of the Boltzmann equation. It is shown that this particular set naturally yields a simple and closed procedure to determine higher order equilibrium distribution functions needed to simulate transcritical flow. The model is validated through several classical benchmarks and is proven to correctly and accurately simulate both 1D and 2D transitions between the two flow regimes. © 2014 Elsevier Inc. All rights reserved.
1. Introduction In the last two decades, the Lattice Boltzmann Method has known a growing popularity for the simulation of a variety of complex flows [3,2]. In particular it has been recognized that some Lattice Boltzmann models can be deduced by discretizing the continuous Boltzmann equation over a velocity space of finite dimension and by defining the Equilibrium Distribution Function (EDF) as a Taylor series expansion of the Maxwellian EDF with respect to the flow velocity [18]. The most commonly employed lattices (e.g. 2DQ9 and 3DQ19) guarantee the equivalence of the lattice Boltzmann model with the Navier–Stokes equation only in the low Mach number limit. This is due to the fact that the velocity sets connected to such lattices are able to exactly reproduce only the lower order hydrodynamic moments of the Maxwellian EDF. In order to correctly reproduce higher order hydrodynamic moments, more complex discrete EDFs are needed which, in turn, require higher dimension velocity spaces [19]. In other words, the highest order of the exactly reproduced hydrodynamic moment is related to the truncation order of the Taylor series expansion of the Maxwellian EDF, which has to be supported by a suitable discretization of the velocity space [27]. A similar limitation is shared by the commonly employed 2DQ9 Lattice Boltzmann models equivalent to the shallow water equations [29]. These models, as for the Navier–Stokes equivalent Lattice Boltzmann models, are derived in the limit of low Froude number, the hydraulic counterpart of the Mach number. It is worth recalling that the Froude number is defined as Fr = U / gh, where U is the characteristic depth-averaged flow velocity, h the water depth, and g the gravitational acceleration. As far as Navier–Stokes equation is concerned, the restriction to low Mach number flows still leaves room for the simulation of a number of technically and theoretically interesting flows. On the contrary, in the shallow water framework, the corresponding limitation on the Froude number is a serious shortcoming for
*
Corresponding author. Tel.: +39 0657333453. E-mail address: [email protected] (M. La Rocca).
http://dx.doi.org/10.1016/j.jcp.2014.12.029 0021-9991/© 2014 Elsevier Inc. All rights reserved.
118
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
real-world applications, where transcritical flows (i.e. flows for which Fr is greater than one) are commonly encountered. A recent work [5] has shown that a 1D transcritical shallow water lattice Boltzmann model can be constructed by means of asymmetrical lattices. The method proposed in [5] has the noteworthy advantage to keep the algorithmic simplicity of the standard lattice Boltzmann models, though permitting the simulation of 1D transcritical shallow water flows. Another way to achieve a Boltzmann-based Finite-Volume numerical model, aimed to solve supercritical shallow water flows, consists in expressing the conservative terms of the shallow water equations as statistical moments of the distribution function with respect to the particle velocity, attaining the so-called Gas Kinetic (hereinafter GK) scheme. This approach has been proposed and validated [8,12,28,9,16], but the intrinsic and intriguing simplicity of the lattice Boltzmann model is lost, due to the rather cumbersome calculations required in order to determine intercell fluxes [8,12]. The aim of this work is to present and validate a generic multispeed 2D shallow water Discrete Boltzmann Model (DBM), which partially recovers the original simplicity of the Lattice Boltzmann Model, while being able to simulate trans- and supercritical 1D and 2D shallow water flows. Hereinafter, by “multispeed”, we mean 2D velocity sets with more than 9 elements. An ideal discretization of the velocity space should be such to allow for an arbitrary number of velocities, comply with the isotropy requirements and preserve exactness of the streaming phase, i.e. generate a space filling lattice. The fact that most multispeed discretizations yield non-space filling lattices [19] is the main shortcoming of such multispeed velocity sets, because the original simplicity of the Lattice Boltzmann method, intrinsically connected to the streaming phase, is lost. In this work this shortcoming is overcome by employing a conventional finite difference scheme for the solution of the resulting multispeed 2D shallow water DBM, although any other discretization method (Finite Volume [15] or Finite Element [11]) could be applied. The structure of the paper is as follows: first, the shallow water equations are recalled for the reader’s convenience; second, it is briefly proved how the classical Gauss–Hermite quadrature cannot be employed to formulate a shallow water Lattice Boltzmann model and thus the derivation of a generic multispeed DBM is shown; third, the ability of the proposed DBM is tested by means of some selected benchmark cases; fourth, results are discussed and conclusions are drawn. 2. Shallow water equations The 2D shallow water equations set is:
∂U +∇ ·E=S ∂t where:
U=
h uh vh
(1)
E=
hu hv hu + gh2 /2 huv huv hv 2 + gh2 /2 2
S=
0 Sx Sy
(2)
and h, u , v are respectively the water depth and depth-averaged velocity components along horizontal directions x and y. In the above, S x and S y are the components of the external force per unit mass along x and y, usually encompassing various effects (bed slope, bed friction, wind-induced surface stress, etc.). The symbol ∇· stands for divergence operator. The homogeneous part of (1) can be obtained from the 3D, free surface formulation of the Euler equations, in the limit of “long waves”, that is, the perturbations of the free surface having a length much larger than the undisturbed depth [13]. 3. Systematic derivation of multispeed EDFs 3.1. Derivation based on Gauss–Hermite quadratures In this section it is shown how the usually employed Gauss–Hermite quadrature cannot be used to derive kinetic models equivalent to the shallow water equations. Consider the following shallow water Maxwellian EDF [8]:
f (h, u, c) =
1
πg
e
− (c−u)·(2 c−u) c
(3)
h
where c = (c x , c y ) is the particle velocity, ch = gh and u = (u , v ). The (n + m)th hydrodynamic moment I nm (n = 0, 1 . . . ; m = 0, 1, . . .) is defined as the statistical moment of (3): +∞ ¨
I nm =
f (h, u, c)cnx cm y dc x dc y
(4)
−∞
I nm is a generic tensor, whose order is n + m. Zero, first and second order hydrodynamic moments I 00 , { I 10 , I 01 }, {{ I 20 , I 11 }, { I 11 , I 02 }} are a scalar, a vector and a second order tensor respectively, namely the water depth, the specific discharge and the momentum flux:
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
119
I 00 = h I 10 = hu I 01 = hv I 20 =
gh2
2 I 11 = huv
I 02 =
gh2 2
+ hu 2
+ hv 2
(5)
Assuming a finite number of particle velocities (hereinafter referred to as velocity set), a set of EDFs is introduced, which can be approximated by a power expansion of the Maxwellian EDF (3) in flow velocity. Standard small particle velocity sets, such as the symmetrical D1Q3 [26] and the D2Q9 [29], derived from a low order (namely second) Taylor expansion of the Maxwellian (3), are not able to simulate transcritical flows. In the following it is shown that, in analogy with what has been shown in literature for Lattice Boltzmann models equivalent to Navier Stokes equations (see [19] for an extensive reference on multispeed models), EDFs derived from high order (namely fourth) Taylor expansion of the Maxwellian (3) allow for the simulation of shallow water transcritical flows. The systematic derivation of such EDFs and the corresponding velocity sets is based on the Taylor series expansion of (3), with respect to u, v:
f (h, u, c) ≈
1
πg ·e
e
−
−
c 2y c2 h
c2 x c2 h
1+2
1+2
cx u ch2
cy v ch2
−
−
u2 ch2
v2 ch2
+2
+2
cx u
2
ch2
cy v
2
ch2
4
+ 4
+
cx u ch2
3
3
3
cy v
−2
3
ch2
−2
u3cx
(ch2 )2
v 3c y
(ch2 )2
+ ...
+ ...
(6)
It is worth noting that the Taylor expansion (6) is performed with respect to the ratio between the fluid velocity u and the linear wave celerity ch = gh, i.e. with respect to the Froude number of the flow. Increasing the Froude number, a corresponding increasing number of terms in the Taylor expansion (6) is needed in order to correctly reproduce the Maxwellian EDF (3). The adoption of an n + m truncation order in (6) permits to exactly reproduce hydrodynamic moments up to the n + m c order. If the scaling ξ = cc x , η = c y is introduced into (6), the expression of the hydrodynamic moment (4) can be factorized h h as follows:
I nm = h
chn+m
π
ˆ ˆ +∞ +∞ n −ξ 2 m −η 2 ξ p (ξ )e dξ η q(η)e dη −∞
(7)
−∞
p, q are real valued functions defined by:
p (ξ ) = 1 + 2 q(η) = 1 + 2
ξu ch
ηv ch
− −
u2 ch2 v2 ch2
+2 +2
2
ξu
+
ch
ηv ch
2 +
4
3 4 3
ξu ch
ηv ch
3 −2
u3 ξ
3 −2
ch3 v 3η ch3
+ ... + ...
(8)
Integrals in (7) can be approximated by the Gauss–Hermite formula [1]:
I nm ≈ hchn+m
N M N ! M ! 2 N + M −2 (ξi )n (η j )m p (ξi )q(η j )
N2 M2
i =1 j =1
( H N −1 (ξi ) H M −1 (η j ))2
(9)
being H N , H M the Nth, Mth order Hermite polynomials and ξi and η j the ith, jth roots of the Nth, Mth order Hermite polynomials respectively. From (9) it is straightforward to obtain the definitions of the velocity vectors and of the corresponding EDFs. Indeed, define the integer k as k = i j. The Cartesian components of the kth velocity vector are obtained as:
ck = ch {ξi , η j } 1 ≤ i ≤ N , 1 ≤ j ≤ M , 1 ≤ k ≤ N M
(10)
As a consequence, (9) takes the form:
I nm ≈ h
N ×M k =0
e cnxk cm yk w k f k
(11)
120
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
where the kth EDF is given by:
f ke = p (ξi )q(η j )
(12)
and the kth weight coefficient, relative to the kth EDF, by:
w k = chn+m
2 N + M −2
N !M !
(13)
( N M )2 ( H N −1 (ξi ) H M −1 (η j ))2
It is interesting to observe that, for any given N and M, the following properties hold for the weights and the velocity vectors (10): N ×M
wk = 1
k =0 N ×M
w k c xk =
k =0 N ×M
w k c yk = 0
k =0
w k c 2xk =
k =0 N ×M
N ×M
N ×M
w k c 2yk =
k =0
ch2 2
w k c xk c yk = 0
(14)
k =0
The third formula in (14) gives the “sound” velocity c s associated to the chosen velocity set: c 2s = ch2 /2 as in Zhou [29]. Eq. (10) allows for the definition of an arbitrary number of velocity sets and of the corresponding EDFs. As stated above, not all such multispeed sets generate Cartesian lattices. But this is not the major shortcoming. Indeed, such a systematic procedure for the definition of sets with an arbitrary number of velocities cannot be adopted in the shallow water framework because the velocity components in (10) depend on the water depth h and thus are not constant. This implies that the set of the Boltzmann–Bhatnagar, Gross and Krook (BGK) kinetic equations:
f e − fk ∂ fk + ck · ∇ f k = k ∂t τ
(1 ≤ k ≤ N × M )
(15)
in force of (10) is not equivalent to the homogeneous part of the shallow water equations (1), being ck dependent on h and thus on space and time. In other words, a space–time dependent sound-speed rules out the use of standard equilibria obtainable by the above Gauss–Hermite quadrature based on the works of Shan et al. [19]. 3.2. Derivation based on the matching of hydrodynamic moments A velocity set with N T + 1 elements and the corresponding set of EDFs, to be used in the shallow water framework, can be defined starting from a polynomial expression of the kth EDF [27]. In the present work, such polynomial expression retains terms up to fourth order:
f ke = h A k + B k
+ Fk
u · ck ch2
+ Ck
u·u ch2
+ Dk
u · ck
2
ch2
+ Ek
u · ck
3
ch2
2 4 (u · u)(u · ck ) u·u (u · u)(u · ck )2 u · ck + G + H + I k k k (ch2 )2 ch2 (ch2 )3 ch2
(16)
It is worth noting that the polynomial form (16) is expressed in terms of the ratio of the fluid velocity u over ch = gh, i.e. of the Froude number of the flow, as the Taylor expansion (6). According to Watari and Tsutahara [27], the unknown coefficients A k , . . . , I k have to be determined by matching the discrete hydrodynamic moments of the EDFs (16) with the exact expressions (4): +∞ ¨ NT − (c−u)·(2 c−u) 1 c m n e n h c xk c yk f k = cm dc x dc y x cye k =0
πg
(17)
−∞
where c xk , c yk are the Cartesian components of the kth velocity ck , k = 0, . . . , N T . It is necessary to obtain a number of matching conditions from (17) equal to the number of unknown coefficients. Depending on the choice of the velocity set, the problem can be underdetermined. However, it is possible to reduce the number of coefficients by exploiting the isotropy properties of the velocity sets [27]. Indeed, velocities having the same magnitude share the same set of coefficients
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
121
Fig. 1. Sketch of the velocity set.
A k , . . . , I k . For this reason, it is useful to group them into subsets, hereinafter referred to as shells, based on their magnitude. In the following we will assume two shells of velocity vectors and the vanishing velocity:
⎧ k=0 ⎨ {0, 0}, 4π 4π ac { cos ( k ), sin ( k )}, 1 ≤ k ≤ N2T 0 ck ≡ NT NT ⎩ bc 0 {cos( 4Nπ k), sin( 4Nπ k)}, N2T < k ≤ N T T T
(18)
being a, b the dimensionless velocity magnitudes of the two shells, scaled by the constant velocity c 0 . In Fig. 1 an example with N T = 16, is shown. The number of unknown coefficients appearing in the definition of the EDFs (16) depends only on the number of shells and on the highest order of terms appearing in (16). Generally speaking, if ns is the number of shells, the definition of EDFs in (16) needs 3 + 9ns unknown coefficients. In Appendix A the expressions of the 21 coefficients relative to the two velocity shells (18) are reported: these expressions, obtained imposing matching conditions (17), are valid for any given a, b, c 0 , N T . It is worth noting that coefficients A k satisfy the same properties (14) of coefficients w k : NT
Ak = 1
k =0 NT
A k c xk =
k =0 NT
A k c yk = 0
k =0
A k c 2xk =
k =0 NT
NT
NT
A k c 2yk =
k =0
ch2 2
A k c xk c yk = 0
(19)
k =0
A Chapman–Enskog expansion, performed according to the work of Sterling and Chen [20], shows that the hydrodynamic limit of the kinetic equations (15) is the set of shallow water equations (2) with an approximation error proportional to 3 , being the smallness parameter of the Chapman–Enskog expansion. It is well known [20] that the latter is the Knudsen number of the flow, expressed as:
=
c0tc
(20)
L
being t c the mesoscopic timescale and L the macroscopic lengthscale. t c has the same order of magnitude of the relaxation time τ , then [22,25]:
tc ∼ τ ∼
ν c 02
being ν the fluid kinematic viscosity. Then the order of magnitude of the smallness parameter the Reynolds and Froude numbers of the flow as follows:
(21)
is expressed in terms of
122
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
∼O
Fr
(22)
Re
Being the Reynolds number of the flow defined as Re = UνL . The calculations of the Chapman–Enskog expansion are standard but rather tedious and are briefly reported in Appendix B for the sake of conciseness. 4. External forces External forces of various types can be introduced as source terms in (15) by means of the following expression:
φk = 2 N
F · ck T
j =0
(23)
( c j · c j)
being F the vector of the external forces. This study only deals with the force induced by the bed slope, whose expression is:
F = − gh∇ zb
(24)
where zb = zb (x, y ) is the bottom elevation. Cartesian components F x , F y of the external force are obtained as:
Fx =
NT
φk c xk
k =0
Fy =
NT
φk c yk
(25)
k =0
Then the kinetic equations (15) with external forces take the form:
f e − fk ∂ fk + ck · ∇ f k = k + φk k = 0, . . . , N T ∂t τ
(26)
5. Results The N T + 1 kinetic equations (26) are solved by means of a conventional finite difference numerical algorithm on a structured staggered 2D uniform Cartesian grid, employing an explicit first order discretization of the time derivative and a first order upwind discretization of the space derivative [10]. It is worth mentioning that any other discretization method, such as Finite Volume [15] or Finite Element [11], can be applied to the N T + 1 kinetic equations (26). As per the treatment of boundary conditions, equilibrium boundary conditions are employed as in Ubertini et al. [24]: the procedure consists of calculating the EDFs with the imposed macroscopic variables and then using them as boundary values for the finite difference scheme. Stability is ensured by keeping the Courant number C = U max t / lower than 1, being U max the highest celerity of perturbations in the fluid domain and = x = y the grid spacing. It has been previously shown that stable Lattice Boltzmann-based shallow water simulations require high values for the relaxation time τ ∗ , and the resulting viscosity can restrict the field of applicability of such models [17]. When solving directly Eq. (26), the relaxation time needs to induce positive viscosity, i.e. τ ∗ > 0 [25]. For the cases considered here the minimum value for the relaxation time τ ∗ ensuring stability always resulted to be lower than one. In order to assess the ability of the multispeed DBM defined by the EDFs (16), the velocity set (18) and the kinetic equations (26) in simulating transcritical and supercritical shallow water flows, the following benchmark cases have been considered: 1) the one-dimensional (1D) dam-break over a flat surface, 2) the 1D steady flow over an uneven bottom profile, 3) the 2D circular dam break over a horizontal bed and 4) the 2D dam break over a horizontal bed. 5.1. 1D dam-break The 1D dam break over a flat bottom is a very simple and yet effective case for assessing accuracy and reliability of any numerical method for the shallow water equations. It deals with the transient evolution of an initial discontinuity of the water level located at x = 0. At t = 0 the water level h is equal to hm , for x ≤ 0, and to h v , for x > 0. This Riemann problem admits the analytical solution of Stoker [21]:
⎧ hm , ⎪ ⎨ 1 (2 ghm − xt )2 , h = 9g∗ ⎪ ⎩h , hv ,
⎧ ⎪ u , x < t ghm ⎪ ⎨ 2m x ∗ ∗ u = 3 ( ghm + t ), t ghm ≤ x < (u − c )t ∗ ∗ ∗ ⎪ (u − c )t ≤ x < st ⎪ ⎩u , uv st ≤ x
(27)
where s is the shock’s propagation speed, which is a function of the intermediate water depth and velocity, h ∗ , u ∗ according to the following relations:
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
123
Fig. 2. 1D dam-break, comparison between analytical (solid traces) and numerical results (markers): upper panel, dimensionless water depth (black traces) and fluid velocity profiles (gray traces); lower panel, Froude number profile. t ∗ = 0.12: ; t ∗ = 0.17: ; t ∗ = 0.29: ; where t ∗ = t ghm / L.
⎧ ⎪ 1+ 1+ ⎪ ⎪ ⎪ ⎨ u ∗ = s − gh v s
8s2 gh v
4s
⎪ h∗ = − ∗ hv ⎪ ⎪ ⎪ u − s ⎩ ∗ ∗ u + 2 gh = 2 ghm
(28)
The typical shape of the Stoker’s solution for the 1D dam-break flow is shown in Fig. 2, where it is compared to the DBM results in terms of depth, flow velocity and Froude number. Numerical results are normalized by hm (water depth), ghm (flow velocity) and L (channel’s length). The agreement is considerably good and it has been quantified by the mean absolute error, whose definition for a generic variable q is:
´ Err[q] =
|qref (x) − q(x)|dx ´ |qref (x)|dx
(29)
For the case considered Err[h] ∼ 10−2 . The agreement between numerical and analytical flow velocity is as good as for the water depth, the corresponding error Err[u ] having the same order of magnitude. The comparison in terms of Froude number shows a discrepancy between numerical and analytical results near the shock which heavily influences the value of a lumped error measure such as (29). The overall agreement is anyway remarkable. It is worth observing the ability of the proposed DBM to smoothly simulate the sub-supercritical transition, without any instability issue, in a case for which the maximum Froude number is equal to: Fr = 5.74. The simulation, though intrinsically 1D, has been performed in a 2D domain discretized by 1000 × 5 nodes. Free-slip boundary conditions are imposed everywhere.
124
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
Three different velocity sets have been implemented, obtained by increasing the number of allowed velocities: 21, 41, 81. It is worth noting that no significant improvement of numerical results has been observed by increasing the number of allowed velocities, as long as it is higher than a minimum value, below which the simulation becomes unstable. It is worth noting that instabilities are generated only if a critical transition (Fr = 1) occurs. The spatial resolution affects the accuracy of the front position, i.e. the lower the height h v the finer the grid needs to be for an accurate representation of the front position, as usual for the numerical solution of shallow water models. Values of the parameters a, b, c 0 have been set to 1, 2, ghm respectively. The minimum value of the relaxation time τ ∗ ensuring stability was 0.8. Courant number was set to 0.1. 5.2. 1D steady flow over a bottom profile In this benchmark the ability of simulating a transcritical flow induced by external forces is tested. More specifically, a steady state solution with a known analytical formulation is considered. Consider a 1D straight channel of length L, described by the x abscissa, and having the following bed elevation profile:
zb (x) = α e −(
x−x1 2
σ
)
+ β e −(
x−x2 2
σ
)
+ γ e −(
x−x3 2
σ
)
(30)
which consists of three consecutive Gaussian bumps, whose crests are respectively located at x1 , x2 , x3 , with elevations α , β , γ , and common variance σ . Fixing the specific (per unit width) discharge q0 over the whole domain, the critical depth can be calculated as:
hc =
q0
√
23 (31)
g
Imposing that the current goes through critical depth over the top of the second and third bumps (i.e. at x = x2 , x = x3 ), it is possible, under the hypothesis of no energy dissipation, to analytically derive the steady water depth profile shown in Fig. 3. This steady motion consists of a subcritical–subcritical flow over the first bump at x = x1 , followed by two subcritical–supercritical transitions over the second and third bumps at x = x2 , x = x3 respectively. Downstream of the first sub-supercritical transition a super-subcritical transition, i.e. a hydraulic jump, occurs. The hydraulic jump occurs where the upstream supercritical specific thrust equals the downstream subcritical one, that is where: F (h u , q0 ) = F (hd , q0 ), the thrust being defined as:
F (h, q0 ) = ρ gh3 + 2q20 /2h
(32)
h u , hd , ρ are the water depths upstream and downstream of the hydraulic jump and the density of water respectively. The steady state DBM numerical profile is obtained evolving from an initial uniform water depth, keeping the upstream water depth and discharge fixed. The numerical solution is considered to be steady when the normalized maximum increment I s between two consecutive timesteps n and (n − 1):
I s = max i
|hni − hni −1 | hni −1
(33)
satisfies the condition: I s < 10−8 . The steady state is reached at t ∼ 100 s. In the upper panel of Fig. 3 the dimensionless analytical and numerical water depth profiles are shown. The water depth scale is the critical depth (31). In the lower panel of Fig. 3 the analytical and numerical Froude number profiles are shown. All profiles in Fig. 3 employ a horizontal length scaling equal to σ . The agreement between numerical and analytical data is remarkably good. For the case considered the error (29) is: Err[h] ∼ 10−3 . The super-subcritical transition occurs as a discontinuity, i.e. as an abrupt elevation, of the water depth. This behavior is reflected in the plot of the Froude number shown in Fig. 3: immediately downstream the smooth sub-supercritical transition occurring over the second bump at x = x2 , the inverse transition is revealed by a sudden decrease of the Froude number, caused by an abrupt increase of the depth and a corresponding decrease of the flow velocity, i.e. a bore or hydraulic jump occurring at 31.3 ≤ x/σ ≤ 31.4. It is worth noting that the bed profile has been chosen in such a way to make this super-subcritical transition occur at a location where the bed slope, though very small, is not null: the analytical solution for both the depth and the flow velocity profile immediately downstream the hydraulic jump thus consists in a smooth variation toward the downstream horizontal profile. The simulation, though intrinsically 1D, has been performed in a 2D domain discretized by 1000 × 5 nodes. Free-slip conditions have been imposed on all boundaries. For this test case, the parameters a, b, c 0 have been set to 1, 2, gh0 respectively. A set of 81 velocities has been employed. The relaxation time was chosen as τ ∗ = 0.8. Courant number was set to 0.1. This test and the previous one demonstrate how the proposed model is able to flawlessly simulate 1D supercritical shallow water regimes with Froude number much higher than one. It is interesting to observe that the numerical steady flow agrees remarkably well with the analytical solution, which has been obtained in absence of dissipative forces. This suggests that the main source of numerical viscosity for the model under consideration stems from the first order finite difference discretization, and not from the viscosity of the underlying kinetic model. This will be shown in detail in Section 6.
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
125
Fig. 3. 1D steady flow profiles over Gaussian bumps, comparison between analytical (solid traces) and numerical (markers) results: upper panel, dimensionless water depth; lower panel, Froude number.
5.3. 2D circular dam-break The 2D circular dam break test of Billett and Toro [4] is here considered. The problem consists in the collapse of a circular cylinder of liquid, of initial depth H and radius R, surrounded by fluid at rest with initial depth h. The radial symmetry can be exploited in order to numerically compute a 1D shallow water reference numerical solution, which depends only on the radial coordinate and time: here a classical Runge Kutta 22 method has been used to integrate the 1D set of equations; the reference solution for the case h/ H = 1/10, R / H = 3.5 was generated discretizing R with 200 computational points, and used to assess DBM accuracy. In Fig. 4 the water depth and Froude number are shown at two different instants of time. Froude number reaches Fr 2 during the simulation, making this test ideal for assessing the capability to simulate 2D transcritical flows. The considered simulation has been run on a 600 × 600 Cartesian grid, Courant number was set to 0.1. From Fig. 4 it is evident that the DBM numerical solution, which has been formulated as 2D, matches remarkably well the reference one. The smearing visible at some locations is imputable to the numerical viscosity and will be discussed in Section 6. Error has been calculated for the Froude number and water depth using a formula similar to (29) but accounting
126
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
Fig. 4. Circular dam break: comparison in terms of Froude number (dashed lines and empty markers, left panel), and water depth (continuous lines and filled markers, right panel), for two instants (gray traces t = 0.2, black traces t = 0.5).
for problem axial symmetry. A value slightly lower than 10−2 for the depth and 10−1 for the Froude number has been found for the case considered. 5.4. 2D dam-break The 2D dam break test of Fennema and Chaudhry [6] is here considered. It consists of the propagation of a wave triggered by the instantaneous collapse of a lock separating two parts of the domain, each one with a different initial water level at rest (see the sketch in Fig. 5 for reference). The square domain has a side of 200 m, is divided into two equal parts by a 10 m thick wall. The wall has a 75 m wide breach, extending from y = 95 m to y = 170 m. The initial water level is equal to hm on the left of the wall (namely x ≤ 95 m), and to h v on the right of the wall (namely x > 95 m). At t = 0 s the breach is considered open and a 2D dam break flow is generated. Numerical simulations described in the following are relative to a case with hm = 10 m, h v = 0.3 m. Such markedly high difference was chosen to produce transcritical 2D shocks. The simulation is characterized by the propagation of a weak, rounded-shaped shock traveling almost perpendicularly to the dividing wall, meanwhile, a rarefaction wave spreads radially into the left part of the domain with a speed which is almost half of the shock’s one ((27) gives a close estimate of such celerities). The shock impacts on the opposite wall and is reflected in the form of a strong, wide and backward traveling hydraulic jump; the accumulation of water along the opposite wall spreads laterally and impinges into both lateral walls creating a complicated system of waves traveling back toward the breach and propagating upwind through the aperture into the still basin. No analytical solution is available for this benchmark thus in this work a well established and validated numerical model has been considered as reference. The reference numerical model integrates the 2D shallow water equations (1) by means of a Finite Volume, shock capturing scheme over an unstructured triangular mesh [14]. It employs a second order Total Variation Diminishing, Weighted Average Flux method [23]. Such scheme has been developed in order to guarantee correct propagation speed of discontinuous solutions, while maintaining a second order accuracy over smooth ones. The same benchmark has been simulated also by means of the GK scheme [16] in order to compare the accuracy and the efficiency of the proposed multispeed approach with another Boltzmann-based algorithm able to simulate supercritical flows. The DBM and the GK simulations have been carried out on a 100 × 100 grid. Values of the DBM parameters a, b, c 0 have been set to 1, 2, ghm respectively. The employed velocity set retained 81 elements and τ ∗ was set to 0.8. Courant number has been set to 0.1. The reference numerical simulation has been carried out on an unstructured mesh of 3574 triangular cells. The assessment of the accuracy of the present model is carried out comparing the time histories of depth and Froude number (Fig. 6), and specific discharges q x , q y (Fig. 7) at points P 1 , P 2 , P 3 , P 4 (see Fig. 5 for points’ locations), with the results yielded by the reference model. Fig. 6 and Fig. 7 include GK results. Water depths are scaled with hm . Specific discharges q x , q y are defined as: q x = uh, q y = vh and scaled with q0 = hm
ghm . Finally, time is scaled with t 0 = hm . ghm
The three models substantially yield the same results, being the shape of the time histories very similar at most of the measuring points. The propagation of steep shocks, observable at P 2 for example, is well reproduced both in terms of strength and speed, with a slight advance of the proposed model for what concerns the bore reflected by the front wall. It is worth noting that the flow becomes markedly supercritical, as can be seen by inspecting the Froude number values attained at points P 1 and P 2 in Fig. 6. The proposed model shows a marked numerical viscosity compared to the reference
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
127
Fig. 5. Sketch of the 2D dam break benchmark of Fennema and Chaudhry [6].
Fig. 6. Time histories of dimensionless water depth (black trace) and Froude number (gray trace) at gauge locations. Comparison among DBM numerical model (solid traces), GK numerical model (dotted traces) and reference solution (dashed traces).
and the GK model, due to both the adopted time–space integration of the governing (26) and the value of the relaxation time. Table 1 shows the values of mean absolute errors Err [h], Err[Fr] (29) between the two models (the DBM and the GK) and the reference one in order to compare the numerical accuracy. The higher accuracy of the proposed DBM with respect to that of the GK model is evident. The main purpose of this benchmark is, as stated above, to check for the ability of the model to correctly compose multidirectional transcritical shocks: based on the outcome of the comparison, it can be concluded that the model possesses such feature. In addition to such quantitative assessment, a general idea on the ability
128
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
Fig. 7. Time histories of dimensionless specific discharge (along x-direction black traces, along y-direction, gray traces) at gauge locations. Comparison among DBM numerical model (solid traces), GK numerical model (dotted traces) and reference solution (dashed traces).
Table 1 Comparison between the mean absolute error attained by DBM and GK numerical results with respect to the reference solution at points P 1 , P 2 , P 3 and P 4 .
P1 P2 P3 P4
DBM–REF Err (h)
GK–REF Err (h)
DBM–REF Err(Fr)
GK–REF Err(Fr)
0.06 0.06 0.03 0.02
0.09 0.06 0.05 0.04
0.11 0.09 0.12 0.26
0.31 0.12 0.26 0.24
Fig. 8. 2D plot of the water depth in the whole domain at t = 26.9 s. Left panel: DBM numerical results. Right panel: numerical reference results. Distances and water depth values are expressed in meters.
of the proposed multispeed DBM can be gleaned from Fig. 8, which shows the distribution of the water surface elevation in the whole domain at t = 26.9 s. In Fig. 8 both the horizontal and the vertical lengthscales are expressed in meters. All the most important flow structures are very similar in both panels, in particular the shape and the position of the curved shock. Finally, as far as the computational performance is concerned, the DBM, due to its intrinsically simple structure, presents a slighter advantage in terms of computational time with respect to both the GK and the reference model. However, further dedicated tests are needed to draw any conclusion on the computational burden needed by the proposed method.
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
129
6. Kinetic and numerical viscosity The Shallow Water Equations are usually expressed as an inviscid set of conservation laws. All the benchmarks employed in this work compare the results of the DBM with inviscid solutions (either analytical or numerical) of the system. The DBM, while showing an overall considerable agreement with those reference solutions, inevitably introduces additional viscosity. The two sources of dissipation are the viscosity term νkin introduced by the collision step in the kinetic equation (26) and the numerical diffusion νFD generated by the Finite Difference scheme employed to solve it. The latter can be expressed as in Fletcher [7]: νFD = xc 0 (1 − C )/2 where C is the Courant number imposed. The ratio between the two viscosities can be then expressed as:
νFD /νkin =
1−C
(34)
2C τ ∗
The expression for the kinetic viscosity is deduced from (21) estimating τ ∼ τ ∗ t. All simulations were run with C = 0.1 and τ ∗ = 0.8, thus resulting in the numerical viscosity being roughly five times greater than the kinetic one. 7. Conclusions In this work a Discrete Boltzmann Model able to solve the 2D transcritical shallow water equations has been developed and validated. The model employs an original discretization of the continuous particle velocity space consisting of two sets of velocities grouped on the basis of their magnitude. The particular structure of the chosen velocity set allows to significantly reduce the number of unknown coefficients of a fourth order polynomial expression of the Equilibrium Distribution Functions: the coefficients are obtained by matching discrete hydrodynamics moments up to fourth order with their continuous counterparts. Several benchmark cases have been carried out in order to allow for a thorough assessment of the ability of the proposed DBM to correctly converge to the solution of the shallow water equations when 1D or 2D supercritical flows develop. This study shows that Boltzmann-based methods can be extended to the simulation of trans- and supercritical shallow water flows, frequently found in real-world applications. Acknowledgement The research has been supported by the Italian national project “Hydroelectric energy by osmosis in coastal areas”, PRIN 2010–2011. Appendix A. Definition of the EDF coefficients
If we define φ = c 0 / gh, for k = 0:
A 0 = 1 − a2 + b 2 −
C0 =
4
φ2
/ φ 2 a2 b 2 2
2
φ
− a2 − b 2 / a2 b 2
G 0 = 1/ a 2 b 2
Bk = β
2(2 − φ 2 b2 ) a2 φ 4 4(2 − φ 2 b2 ) a2 φ 2
C k = −β
Dk = β
Ek = β
2(2 − φ 2 b2 ) a2 φ 2
8(3 − φ 2 b2 ) a4 φ 2 16 3a4
(A.2)
Defining β = 1/[ N T (a2 − b2 )], for 1 ≤ k ≤
Ak = β
(A.1)
(A.3) NT 2
:
(A.4)
(A.5)
(A.6)
(A.7)
(A.8)
130
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
Fk = 0
(A.9)
G k = −β
1
(A.10)
a2
Hk = 0 Ik = β For
NT 2
(A.11) 8
(A.12)
a6
+ 1 ≤ k ≤ NT : 2(a2 φ 2 − 2)
Ak = β
4(a2 φ 2 − 2)
Bk = β
2(a2 φ 2 − 2) b2 φ 2
8(a2 φ 2 − 3)
Dk = β
16(2a2 − 3b2 ) 3b6 8(a2 − b2 )
b4
H k = −β
24(a2 − b2 )
8 b8
(A.17)
(A.18)
b4
1 2 3a − 2b2
Gk = β
(A.15)
(A.16)
b4 φ 2
F k = −β
Ik = β
(A.14)
b2 φ 2
C k = −β
Ek = β
(A.13)
b2 φ 4
b6 3a2 − 4b2
(A.19)
(A.20)
(A.21)
Appendix B. Chapman–Enskog expansion of the kinetic equations (26) Consider the scaling:
t=
c0t L
= L∇ ∇ τ=
= ck =
τ
tc c0tc L
ck c0
fk fk = , H
fe = k
f ke H
(B.1)
where H , L, t c are respectively the macroscopic vertical and horizontal lengthscales and the mesoscopic timescale, i.e. the time interval between two successive collisions. Applying the scaling (B.1) to time and space derivatives in (26) and expanding them in perturbative series with respect to the smallness parameter (namely the Knudsen number, provided a convective scaling is assumed), the following expressions can be obtained:
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
131
fk = f e + f 1 + 2 f 2 + 3 f 3 + 4 f 4 + . . . k
k
k
k
k
∂ ∂ ∂ ∂ ∂ = + + 2 + 3 + ... ∂ t ∂ t ∂ t ∂ t4 ∂t 1 2 3 = ∇1 + ∇2 + 2 ∇3 + 3 ∇4 + . . . ∇
(B.2)
The following equations are obtained at each order:
e ∂ f ke f1 + ∇1 · ck fk = − k ∂ t1 τ e e ∂ ∂ f k1 fk f2 =− k + ∇1 · ck f k1 + + ∇2 · ck fk ∂ t1 ∂ t2 τ 1 e 2 e ∂ fk f3 ∂ fk ∂ fk =− k + ∇1 · ck f k2 + + ∇2 · ck f k1 + + ∇3 · ck fk ∂ t1 ∂ t2 ∂ t3 τ 2 1 e 3 e ∂ fk f4 ∂ fk ∂ fk ∂ fk 3 2 1 + + + ∇1 · ck f k + + ∇2 · ck f k + ∇3 · ck f k + ∇4 · ck f k = − k ∂ t1 ∂ t2 ∂ t3 ∂ t4 τ
0) 1) 2) 3) .. .
(B.3)
Summing (B.3) with respect to k and accounting for the fact that: NT
h fe = k
k =0
H
NT hu ck f ke =
c0 H
k =0 NT
f ki = 0;
i = 1, 2, . . .
k =0 NT ck f ki = 0;
i = 1, 2, . . .
(B.4)
k =0
it is straightforward to obtain the shallow water mass balance equation, which in dimensional form is given by:
∂h + ∇(hu) = 0 ∂t
(B.5)
Multiplying each (B.3) by ck , summing with respect to k and accounting for the definitions of the hydrodynamic moments:
2 =
NT h2 1 ck ck f ke = 2 hu ⊗ u + g I k =0
3 =
c0 H
2
NT ck ck ck f ke k =0
NT 4 = ck ck ck ck f ke k =0
where symbol ⊗ indicates the dyadic product of the vector u with itself and I is the identity matrix, the following form of the macroscopic momentum balance equation is obtained:
c 02 H h2 ∂ hu ∂ 2 + ∇ · hu ⊗ u + g I = τ ∇1 · + ∇1 · 3 ∂t 2 L ∂ t1 ∂ 2 ∂ 2 2 2 + ∇1 · 3 + τ ∇1 · + ∇2 · 3 + τ ∇2 · ∂ t1 ∂ t2 2 ∂ ∂ 2 3 + τ ∇1 · + 2 ∇ · + ∇ · ∇ · . . . − 2 1 3 1 1 4 ∂ t1 ∂ t 12
(B.6)
132
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
Thus (26), together with the velocity set (18) and the EDFs (16) are equivalent to the shallow water momentum balance equation (2). The equivalence is meant in the limit of small . The approximation error is proportional to 3 , due to the fact that the proposed EDFs (16) together with the velocity set (18) are able to reproduce exactly all the hydrodynamic moments appearing at right hand side of (B.6): i.e. the first term which is not correctly reproduced in (B.6) is proportional to 3 . References [1] Milton Abramowitz, Irene A. Stegun, Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, Courier Dover Publications, 2012. [2] Cyrus K. Aidun, Jonathan R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (1) (2010) 439–472. [3] R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: theory and applications, Phys. Rep. 222 (3) (1992) 145–197. [4] S.J. Billett, E.F. Toro, On WAF-type schemes for multidimensional hyperbolic conservation laws, J. Comput. Phys. 130 (1) (1997) 1–24. [5] Bastien Chopard, V.T. Pham, Laurent Lefevre, Asymmetric lattice Boltzmann model for shallow water flows, Comput. Fluids 88 (2013) 225–231. [6] Robert J. Fennema, M. Hanif Chaudhry, Explicit methods for 2-D transient free surface flows, J. Hydraul. Eng. 116 (8) (1990) 1013–1034. [7] C.A.J. Fletcher, Computational Techniques for Fluid Dynamics, Springer, Berlin, 1991. [8] M.S. Ghidaoui, J.Q. Deng, W.G. Gray, K. Xu, A Boltzmann based model for open channel flows, Int. J. Numer. Methods Fluids 35 (4) (2001) 449–494. [9] M.S. Ghidaoui, A.A. Kolyshkin, J.H. Liang, F.C. Chan, Q. Li, K. Xu, Linear and nonlinear analysis of shallow wakes, J. Fluid Mech. 548 (2006) 309–340. [10] Charles Hirsch, Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics, vol. 1, Butterworth– Heinemann, 2007. [11] Taehun Lee, Ching-Long Lin, A characteristic Galerkin method for discrete Boltzmann equation, J. Comput. Phys. 171 (1) (2001) 336–356. [12] J.H. Liang, M.S. Ghidaoui, J.Q. Deng, W.G. Gray, A Boltzmann-based finite volume algorithm for surface water flows on cells of arbitrary shapes, J. Hydraul. Res. 45 (2) (2007) 147–164. [13] Chiang C. Mei, The Applied Dynamics of Ocean Surface Waves, vol. 1, World Scientific, 1989. [14] Pietro Prestininzi, Numerical modelling of fluvial floods, PhD thesis, Università degli Studi Roma TRE, 2009. [15] Pietro Prestininzi, Michele La Rocca, Reinhard Hinkelmann, et al., Comparative study of a Boltzmann-based finite volume and a lattice Boltzmann model for shallow water flows in complex domains, Int. J. Offshore Polar Eng. 24 (03) (2014) 161–167. [16] Pietro Prestininzi, Michele La Rocca, Andrea Montessori, Giampiero Sciortino, A gas-kinetic model for 2D transcritical shallow water flows propagating over dry bed, Comput. Math. Appl. 68 (4) (2014) 439–453. [17] Pietro Prestininzi, Giampiero Sciortino, Michele La Rocca, On the effect of the intrinsic viscosity in a two-layer shallow water lattice Boltzmann model of axisymmetric density currents, J. Hydraul. Res. 51 (6) (2013) 668–680. [18] X.W. Shan, X.Y. He, Discretization of the velocity space in the solution of the Boltzmann equation, Phys. Rev. Lett. 80 (1) (1998) 65–68. [19] Xiaowen Shan, Xue-Feng Yuan, Hudong Chen, Kinetic theory representation of hydrodynamics: a way beyond the Navier–Stokes equation, J. Fluid Mech. 550 (2006) 413–441. [20] James D. Sterling, Shiyi Chen, Stability analysis of lattice Boltzmann methods, J. Comput. Phys. 123 (1) (1996) 196–206. [21] James Johnston Stoker, Water Waves: The Mathematical Theory with Applications, vol. 36, John Wiley & Sons, 2011. [22] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford Science Publications, Oxford, 2001. [23] E.F. Toro, A weighted average flux method for hyperbolic conservation laws, Proc. R. Soc. Lond. A 423 (1865) (1989) 401–418. [24] S. Ubertini, G. Bella, S. Succi, Lattice Boltzmann method on unstructured grids: further developments, Phys. Rev. E 68 (1) (2003) 1–10. [25] Stefano Ubertini, G. Bella, Sauro Succi, Unstructured lattice Boltzmann equation with memory, Math. Comput. Simul. 72 (2) (2006) 237–241. [26] Pham Van Thang, Bastien Chopard, Laurent Lefèvre, Diemer Anda Ondo, Eduardo Mendes, Study of the 1D lattice Boltzmann shallow water equation and its coupling to build a canal network, J. Comput. Phys. 229 (19) (2010) 7373–7400. [27] Minoru Watari, Michihisa Tsutahara, Possibility of constructing a multispeed Bhatnagar–Gross–Krook thermal model of the lattice Boltzmann method, Phys. Rev. E 70 (1) (2004) 016703. [28] Kun Xu, A well-balanced gas-kinetic scheme for the shallow-water equations with source terms, J. Comput. Phys. 178 (2) (2002) 533–562. [29] Jian Zhou, Lattice Boltzmann Methods for Shallow Water Flows, Springer, Berlin, 2004.
Contents lists available at ScienceDirect
Journal of Computational Physics www.elsevier.com/locate/jcp
A multispeed Discrete Boltzmann Model for transcritical 2D shallow water flows Michele La Rocca a,∗ , Andrea Montessori a , Pietro Prestininzi a , Sauro Succi b a b
Dipartimento di Ingegneria, Universita’ degli Studi Roma TRE, Via Vito Volterra 62, 00146 Rome, Italy Istituto per le Applicazioni del Calcolo, CNR, Via dei Taurini 19, 00189 Rome, Italy
a r t i c l e
i n f o
Article history: Received 23 June 2014 Received in revised form 29 October 2014 Accepted 16 December 2014 Available online 23 December 2014 Keywords: Shallow water equations Multispeed Discrete Boltzmann Model Transcritical flows
a b s t r a c t In this work a Discrete Boltzmann Model for the solution of transcritical 2D shallow water flows is presented and validated. In order to provide the model with transcritical capabilities, a particular multispeed velocity set has been employed for the discretization of the Boltzmann equation. It is shown that this particular set naturally yields a simple and closed procedure to determine higher order equilibrium distribution functions needed to simulate transcritical flow. The model is validated through several classical benchmarks and is proven to correctly and accurately simulate both 1D and 2D transitions between the two flow regimes. © 2014 Elsevier Inc. All rights reserved.
1. Introduction In the last two decades, the Lattice Boltzmann Method has known a growing popularity for the simulation of a variety of complex flows [3,2]. In particular it has been recognized that some Lattice Boltzmann models can be deduced by discretizing the continuous Boltzmann equation over a velocity space of finite dimension and by defining the Equilibrium Distribution Function (EDF) as a Taylor series expansion of the Maxwellian EDF with respect to the flow velocity [18]. The most commonly employed lattices (e.g. 2DQ9 and 3DQ19) guarantee the equivalence of the lattice Boltzmann model with the Navier–Stokes equation only in the low Mach number limit. This is due to the fact that the velocity sets connected to such lattices are able to exactly reproduce only the lower order hydrodynamic moments of the Maxwellian EDF. In order to correctly reproduce higher order hydrodynamic moments, more complex discrete EDFs are needed which, in turn, require higher dimension velocity spaces [19]. In other words, the highest order of the exactly reproduced hydrodynamic moment is related to the truncation order of the Taylor series expansion of the Maxwellian EDF, which has to be supported by a suitable discretization of the velocity space [27]. A similar limitation is shared by the commonly employed 2DQ9 Lattice Boltzmann models equivalent to the shallow water equations [29]. These models, as for the Navier–Stokes equivalent Lattice Boltzmann models, are derived in the limit of low Froude number, the hydraulic counterpart of the Mach number. It is worth recalling that the Froude number is defined as Fr = U / gh, where U is the characteristic depth-averaged flow velocity, h the water depth, and g the gravitational acceleration. As far as Navier–Stokes equation is concerned, the restriction to low Mach number flows still leaves room for the simulation of a number of technically and theoretically interesting flows. On the contrary, in the shallow water framework, the corresponding limitation on the Froude number is a serious shortcoming for
*
Corresponding author. Tel.: +39 0657333453. E-mail address: [email protected] (M. La Rocca).
http://dx.doi.org/10.1016/j.jcp.2014.12.029 0021-9991/© 2014 Elsevier Inc. All rights reserved.
118
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
real-world applications, where transcritical flows (i.e. flows for which Fr is greater than one) are commonly encountered. A recent work [5] has shown that a 1D transcritical shallow water lattice Boltzmann model can be constructed by means of asymmetrical lattices. The method proposed in [5] has the noteworthy advantage to keep the algorithmic simplicity of the standard lattice Boltzmann models, though permitting the simulation of 1D transcritical shallow water flows. Another way to achieve a Boltzmann-based Finite-Volume numerical model, aimed to solve supercritical shallow water flows, consists in expressing the conservative terms of the shallow water equations as statistical moments of the distribution function with respect to the particle velocity, attaining the so-called Gas Kinetic (hereinafter GK) scheme. This approach has been proposed and validated [8,12,28,9,16], but the intrinsic and intriguing simplicity of the lattice Boltzmann model is lost, due to the rather cumbersome calculations required in order to determine intercell fluxes [8,12]. The aim of this work is to present and validate a generic multispeed 2D shallow water Discrete Boltzmann Model (DBM), which partially recovers the original simplicity of the Lattice Boltzmann Model, while being able to simulate trans- and supercritical 1D and 2D shallow water flows. Hereinafter, by “multispeed”, we mean 2D velocity sets with more than 9 elements. An ideal discretization of the velocity space should be such to allow for an arbitrary number of velocities, comply with the isotropy requirements and preserve exactness of the streaming phase, i.e. generate a space filling lattice. The fact that most multispeed discretizations yield non-space filling lattices [19] is the main shortcoming of such multispeed velocity sets, because the original simplicity of the Lattice Boltzmann method, intrinsically connected to the streaming phase, is lost. In this work this shortcoming is overcome by employing a conventional finite difference scheme for the solution of the resulting multispeed 2D shallow water DBM, although any other discretization method (Finite Volume [15] or Finite Element [11]) could be applied. The structure of the paper is as follows: first, the shallow water equations are recalled for the reader’s convenience; second, it is briefly proved how the classical Gauss–Hermite quadrature cannot be employed to formulate a shallow water Lattice Boltzmann model and thus the derivation of a generic multispeed DBM is shown; third, the ability of the proposed DBM is tested by means of some selected benchmark cases; fourth, results are discussed and conclusions are drawn. 2. Shallow water equations The 2D shallow water equations set is:
∂U +∇ ·E=S ∂t where:
U=
h uh vh
(1)
E=
hu hv hu + gh2 /2 huv huv hv 2 + gh2 /2 2
S=
0 Sx Sy
(2)
and h, u , v are respectively the water depth and depth-averaged velocity components along horizontal directions x and y. In the above, S x and S y are the components of the external force per unit mass along x and y, usually encompassing various effects (bed slope, bed friction, wind-induced surface stress, etc.). The symbol ∇· stands for divergence operator. The homogeneous part of (1) can be obtained from the 3D, free surface formulation of the Euler equations, in the limit of “long waves”, that is, the perturbations of the free surface having a length much larger than the undisturbed depth [13]. 3. Systematic derivation of multispeed EDFs 3.1. Derivation based on Gauss–Hermite quadratures In this section it is shown how the usually employed Gauss–Hermite quadrature cannot be used to derive kinetic models equivalent to the shallow water equations. Consider the following shallow water Maxwellian EDF [8]:
f (h, u, c) =
1
πg
e
− (c−u)·(2 c−u) c
(3)
h
where c = (c x , c y ) is the particle velocity, ch = gh and u = (u , v ). The (n + m)th hydrodynamic moment I nm (n = 0, 1 . . . ; m = 0, 1, . . .) is defined as the statistical moment of (3): +∞ ¨
I nm =
f (h, u, c)cnx cm y dc x dc y
(4)
−∞
I nm is a generic tensor, whose order is n + m. Zero, first and second order hydrodynamic moments I 00 , { I 10 , I 01 }, {{ I 20 , I 11 }, { I 11 , I 02 }} are a scalar, a vector and a second order tensor respectively, namely the water depth, the specific discharge and the momentum flux:
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
119
I 00 = h I 10 = hu I 01 = hv I 20 =
gh2
2 I 11 = huv
I 02 =
gh2 2
+ hu 2
+ hv 2
(5)
Assuming a finite number of particle velocities (hereinafter referred to as velocity set), a set of EDFs is introduced, which can be approximated by a power expansion of the Maxwellian EDF (3) in flow velocity. Standard small particle velocity sets, such as the symmetrical D1Q3 [26] and the D2Q9 [29], derived from a low order (namely second) Taylor expansion of the Maxwellian (3), are not able to simulate transcritical flows. In the following it is shown that, in analogy with what has been shown in literature for Lattice Boltzmann models equivalent to Navier Stokes equations (see [19] for an extensive reference on multispeed models), EDFs derived from high order (namely fourth) Taylor expansion of the Maxwellian (3) allow for the simulation of shallow water transcritical flows. The systematic derivation of such EDFs and the corresponding velocity sets is based on the Taylor series expansion of (3), with respect to u, v:
f (h, u, c) ≈
1
πg ·e
e
−
−
c 2y c2 h
c2 x c2 h
1+2
1+2
cx u ch2
cy v ch2
−
−
u2 ch2
v2 ch2
+2
+2
cx u
2
ch2
cy v
2
ch2
4
+ 4
+
cx u ch2
3
3
3
cy v
−2
3
ch2
−2
u3cx
(ch2 )2
v 3c y
(ch2 )2
+ ...
+ ...
(6)
It is worth noting that the Taylor expansion (6) is performed with respect to the ratio between the fluid velocity u and the linear wave celerity ch = gh, i.e. with respect to the Froude number of the flow. Increasing the Froude number, a corresponding increasing number of terms in the Taylor expansion (6) is needed in order to correctly reproduce the Maxwellian EDF (3). The adoption of an n + m truncation order in (6) permits to exactly reproduce hydrodynamic moments up to the n + m c order. If the scaling ξ = cc x , η = c y is introduced into (6), the expression of the hydrodynamic moment (4) can be factorized h h as follows:
I nm = h
chn+m
π
ˆ ˆ +∞ +∞ n −ξ 2 m −η 2 ξ p (ξ )e dξ η q(η)e dη −∞
(7)
−∞
p, q are real valued functions defined by:
p (ξ ) = 1 + 2 q(η) = 1 + 2
ξu ch
ηv ch
− −
u2 ch2 v2 ch2
+2 +2
2
ξu
+
ch
ηv ch
2 +
4
3 4 3
ξu ch
ηv ch
3 −2
u3 ξ
3 −2
ch3 v 3η ch3
+ ... + ...
(8)
Integrals in (7) can be approximated by the Gauss–Hermite formula [1]:
I nm ≈ hchn+m
N M N ! M ! 2 N + M −2 (ξi )n (η j )m p (ξi )q(η j )
N2 M2
i =1 j =1
( H N −1 (ξi ) H M −1 (η j ))2
(9)
being H N , H M the Nth, Mth order Hermite polynomials and ξi and η j the ith, jth roots of the Nth, Mth order Hermite polynomials respectively. From (9) it is straightforward to obtain the definitions of the velocity vectors and of the corresponding EDFs. Indeed, define the integer k as k = i j. The Cartesian components of the kth velocity vector are obtained as:
ck = ch {ξi , η j } 1 ≤ i ≤ N , 1 ≤ j ≤ M , 1 ≤ k ≤ N M
(10)
As a consequence, (9) takes the form:
I nm ≈ h
N ×M k =0
e cnxk cm yk w k f k
(11)
120
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
where the kth EDF is given by:
f ke = p (ξi )q(η j )
(12)
and the kth weight coefficient, relative to the kth EDF, by:
w k = chn+m
2 N + M −2
N !M !
(13)
( N M )2 ( H N −1 (ξi ) H M −1 (η j ))2
It is interesting to observe that, for any given N and M, the following properties hold for the weights and the velocity vectors (10): N ×M
wk = 1
k =0 N ×M
w k c xk =
k =0 N ×M
w k c yk = 0
k =0
w k c 2xk =
k =0 N ×M
N ×M
N ×M
w k c 2yk =
k =0
ch2 2
w k c xk c yk = 0
(14)
k =0
The third formula in (14) gives the “sound” velocity c s associated to the chosen velocity set: c 2s = ch2 /2 as in Zhou [29]. Eq. (10) allows for the definition of an arbitrary number of velocity sets and of the corresponding EDFs. As stated above, not all such multispeed sets generate Cartesian lattices. But this is not the major shortcoming. Indeed, such a systematic procedure for the definition of sets with an arbitrary number of velocities cannot be adopted in the shallow water framework because the velocity components in (10) depend on the water depth h and thus are not constant. This implies that the set of the Boltzmann–Bhatnagar, Gross and Krook (BGK) kinetic equations:
f e − fk ∂ fk + ck · ∇ f k = k ∂t τ
(1 ≤ k ≤ N × M )
(15)
in force of (10) is not equivalent to the homogeneous part of the shallow water equations (1), being ck dependent on h and thus on space and time. In other words, a space–time dependent sound-speed rules out the use of standard equilibria obtainable by the above Gauss–Hermite quadrature based on the works of Shan et al. [19]. 3.2. Derivation based on the matching of hydrodynamic moments A velocity set with N T + 1 elements and the corresponding set of EDFs, to be used in the shallow water framework, can be defined starting from a polynomial expression of the kth EDF [27]. In the present work, such polynomial expression retains terms up to fourth order:
f ke = h A k + B k
+ Fk
u · ck ch2
+ Ck
u·u ch2
+ Dk
u · ck
2
ch2
+ Ek
u · ck
3
ch2
2 4 (u · u)(u · ck ) u·u (u · u)(u · ck )2 u · ck + G + H + I k k k (ch2 )2 ch2 (ch2 )3 ch2
(16)
It is worth noting that the polynomial form (16) is expressed in terms of the ratio of the fluid velocity u over ch = gh, i.e. of the Froude number of the flow, as the Taylor expansion (6). According to Watari and Tsutahara [27], the unknown coefficients A k , . . . , I k have to be determined by matching the discrete hydrodynamic moments of the EDFs (16) with the exact expressions (4): +∞ ¨ NT − (c−u)·(2 c−u) 1 c m n e n h c xk c yk f k = cm dc x dc y x cye k =0
πg
(17)
−∞
where c xk , c yk are the Cartesian components of the kth velocity ck , k = 0, . . . , N T . It is necessary to obtain a number of matching conditions from (17) equal to the number of unknown coefficients. Depending on the choice of the velocity set, the problem can be underdetermined. However, it is possible to reduce the number of coefficients by exploiting the isotropy properties of the velocity sets [27]. Indeed, velocities having the same magnitude share the same set of coefficients
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
121
Fig. 1. Sketch of the velocity set.
A k , . . . , I k . For this reason, it is useful to group them into subsets, hereinafter referred to as shells, based on their magnitude. In the following we will assume two shells of velocity vectors and the vanishing velocity:
⎧ k=0 ⎨ {0, 0}, 4π 4π ac { cos ( k ), sin ( k )}, 1 ≤ k ≤ N2T 0 ck ≡ NT NT ⎩ bc 0 {cos( 4Nπ k), sin( 4Nπ k)}, N2T < k ≤ N T T T
(18)
being a, b the dimensionless velocity magnitudes of the two shells, scaled by the constant velocity c 0 . In Fig. 1 an example with N T = 16, is shown. The number of unknown coefficients appearing in the definition of the EDFs (16) depends only on the number of shells and on the highest order of terms appearing in (16). Generally speaking, if ns is the number of shells, the definition of EDFs in (16) needs 3 + 9ns unknown coefficients. In Appendix A the expressions of the 21 coefficients relative to the two velocity shells (18) are reported: these expressions, obtained imposing matching conditions (17), are valid for any given a, b, c 0 , N T . It is worth noting that coefficients A k satisfy the same properties (14) of coefficients w k : NT
Ak = 1
k =0 NT
A k c xk =
k =0 NT
A k c yk = 0
k =0
A k c 2xk =
k =0 NT
NT
NT
A k c 2yk =
k =0
ch2 2
A k c xk c yk = 0
(19)
k =0
A Chapman–Enskog expansion, performed according to the work of Sterling and Chen [20], shows that the hydrodynamic limit of the kinetic equations (15) is the set of shallow water equations (2) with an approximation error proportional to 3 , being the smallness parameter of the Chapman–Enskog expansion. It is well known [20] that the latter is the Knudsen number of the flow, expressed as:
=
c0tc
(20)
L
being t c the mesoscopic timescale and L the macroscopic lengthscale. t c has the same order of magnitude of the relaxation time τ , then [22,25]:
tc ∼ τ ∼
ν c 02
being ν the fluid kinematic viscosity. Then the order of magnitude of the smallness parameter the Reynolds and Froude numbers of the flow as follows:
(21)
is expressed in terms of
122
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
∼O
Fr
(22)
Re
Being the Reynolds number of the flow defined as Re = UνL . The calculations of the Chapman–Enskog expansion are standard but rather tedious and are briefly reported in Appendix B for the sake of conciseness. 4. External forces External forces of various types can be introduced as source terms in (15) by means of the following expression:
φk = 2 N
F · ck T
j =0
(23)
( c j · c j)
being F the vector of the external forces. This study only deals with the force induced by the bed slope, whose expression is:
F = − gh∇ zb
(24)
where zb = zb (x, y ) is the bottom elevation. Cartesian components F x , F y of the external force are obtained as:
Fx =
NT
φk c xk
k =0
Fy =
NT
φk c yk
(25)
k =0
Then the kinetic equations (15) with external forces take the form:
f e − fk ∂ fk + ck · ∇ f k = k + φk k = 0, . . . , N T ∂t τ
(26)
5. Results The N T + 1 kinetic equations (26) are solved by means of a conventional finite difference numerical algorithm on a structured staggered 2D uniform Cartesian grid, employing an explicit first order discretization of the time derivative and a first order upwind discretization of the space derivative [10]. It is worth mentioning that any other discretization method, such as Finite Volume [15] or Finite Element [11], can be applied to the N T + 1 kinetic equations (26). As per the treatment of boundary conditions, equilibrium boundary conditions are employed as in Ubertini et al. [24]: the procedure consists of calculating the EDFs with the imposed macroscopic variables and then using them as boundary values for the finite difference scheme. Stability is ensured by keeping the Courant number C = U max t / lower than 1, being U max the highest celerity of perturbations in the fluid domain and = x = y the grid spacing. It has been previously shown that stable Lattice Boltzmann-based shallow water simulations require high values for the relaxation time τ ∗ , and the resulting viscosity can restrict the field of applicability of such models [17]. When solving directly Eq. (26), the relaxation time needs to induce positive viscosity, i.e. τ ∗ > 0 [25]. For the cases considered here the minimum value for the relaxation time τ ∗ ensuring stability always resulted to be lower than one. In order to assess the ability of the multispeed DBM defined by the EDFs (16), the velocity set (18) and the kinetic equations (26) in simulating transcritical and supercritical shallow water flows, the following benchmark cases have been considered: 1) the one-dimensional (1D) dam-break over a flat surface, 2) the 1D steady flow over an uneven bottom profile, 3) the 2D circular dam break over a horizontal bed and 4) the 2D dam break over a horizontal bed. 5.1. 1D dam-break The 1D dam break over a flat bottom is a very simple and yet effective case for assessing accuracy and reliability of any numerical method for the shallow water equations. It deals with the transient evolution of an initial discontinuity of the water level located at x = 0. At t = 0 the water level h is equal to hm , for x ≤ 0, and to h v , for x > 0. This Riemann problem admits the analytical solution of Stoker [21]:
⎧ hm , ⎪ ⎨ 1 (2 ghm − xt )2 , h = 9g∗ ⎪ ⎩h , hv ,
⎧ ⎪ u , x < t ghm ⎪ ⎨ 2m x ∗ ∗ u = 3 ( ghm + t ), t ghm ≤ x < (u − c )t ∗ ∗ ∗ ⎪ (u − c )t ≤ x < st ⎪ ⎩u , uv st ≤ x
(27)
where s is the shock’s propagation speed, which is a function of the intermediate water depth and velocity, h ∗ , u ∗ according to the following relations:
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
123
Fig. 2. 1D dam-break, comparison between analytical (solid traces) and numerical results (markers): upper panel, dimensionless water depth (black traces) and fluid velocity profiles (gray traces); lower panel, Froude number profile. t ∗ = 0.12: ; t ∗ = 0.17: ; t ∗ = 0.29: ; where t ∗ = t ghm / L.
⎧ ⎪ 1+ 1+ ⎪ ⎪ ⎪ ⎨ u ∗ = s − gh v s
8s2 gh v
4s
⎪ h∗ = − ∗ hv ⎪ ⎪ ⎪ u − s ⎩ ∗ ∗ u + 2 gh = 2 ghm
(28)
The typical shape of the Stoker’s solution for the 1D dam-break flow is shown in Fig. 2, where it is compared to the DBM results in terms of depth, flow velocity and Froude number. Numerical results are normalized by hm (water depth), ghm (flow velocity) and L (channel’s length). The agreement is considerably good and it has been quantified by the mean absolute error, whose definition for a generic variable q is:
´ Err[q] =
|qref (x) − q(x)|dx ´ |qref (x)|dx
(29)
For the case considered Err[h] ∼ 10−2 . The agreement between numerical and analytical flow velocity is as good as for the water depth, the corresponding error Err[u ] having the same order of magnitude. The comparison in terms of Froude number shows a discrepancy between numerical and analytical results near the shock which heavily influences the value of a lumped error measure such as (29). The overall agreement is anyway remarkable. It is worth observing the ability of the proposed DBM to smoothly simulate the sub-supercritical transition, without any instability issue, in a case for which the maximum Froude number is equal to: Fr = 5.74. The simulation, though intrinsically 1D, has been performed in a 2D domain discretized by 1000 × 5 nodes. Free-slip boundary conditions are imposed everywhere.
124
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
Three different velocity sets have been implemented, obtained by increasing the number of allowed velocities: 21, 41, 81. It is worth noting that no significant improvement of numerical results has been observed by increasing the number of allowed velocities, as long as it is higher than a minimum value, below which the simulation becomes unstable. It is worth noting that instabilities are generated only if a critical transition (Fr = 1) occurs. The spatial resolution affects the accuracy of the front position, i.e. the lower the height h v the finer the grid needs to be for an accurate representation of the front position, as usual for the numerical solution of shallow water models. Values of the parameters a, b, c 0 have been set to 1, 2, ghm respectively. The minimum value of the relaxation time τ ∗ ensuring stability was 0.8. Courant number was set to 0.1. 5.2. 1D steady flow over a bottom profile In this benchmark the ability of simulating a transcritical flow induced by external forces is tested. More specifically, a steady state solution with a known analytical formulation is considered. Consider a 1D straight channel of length L, described by the x abscissa, and having the following bed elevation profile:
zb (x) = α e −(
x−x1 2
σ
)
+ β e −(
x−x2 2
σ
)
+ γ e −(
x−x3 2
σ
)
(30)
which consists of three consecutive Gaussian bumps, whose crests are respectively located at x1 , x2 , x3 , with elevations α , β , γ , and common variance σ . Fixing the specific (per unit width) discharge q0 over the whole domain, the critical depth can be calculated as:
hc =
q0
√
23 (31)
g
Imposing that the current goes through critical depth over the top of the second and third bumps (i.e. at x = x2 , x = x3 ), it is possible, under the hypothesis of no energy dissipation, to analytically derive the steady water depth profile shown in Fig. 3. This steady motion consists of a subcritical–subcritical flow over the first bump at x = x1 , followed by two subcritical–supercritical transitions over the second and third bumps at x = x2 , x = x3 respectively. Downstream of the first sub-supercritical transition a super-subcritical transition, i.e. a hydraulic jump, occurs. The hydraulic jump occurs where the upstream supercritical specific thrust equals the downstream subcritical one, that is where: F (h u , q0 ) = F (hd , q0 ), the thrust being defined as:
F (h, q0 ) = ρ gh3 + 2q20 /2h
(32)
h u , hd , ρ are the water depths upstream and downstream of the hydraulic jump and the density of water respectively. The steady state DBM numerical profile is obtained evolving from an initial uniform water depth, keeping the upstream water depth and discharge fixed. The numerical solution is considered to be steady when the normalized maximum increment I s between two consecutive timesteps n and (n − 1):
I s = max i
|hni − hni −1 | hni −1
(33)
satisfies the condition: I s < 10−8 . The steady state is reached at t ∼ 100 s. In the upper panel of Fig. 3 the dimensionless analytical and numerical water depth profiles are shown. The water depth scale is the critical depth (31). In the lower panel of Fig. 3 the analytical and numerical Froude number profiles are shown. All profiles in Fig. 3 employ a horizontal length scaling equal to σ . The agreement between numerical and analytical data is remarkably good. For the case considered the error (29) is: Err[h] ∼ 10−3 . The super-subcritical transition occurs as a discontinuity, i.e. as an abrupt elevation, of the water depth. This behavior is reflected in the plot of the Froude number shown in Fig. 3: immediately downstream the smooth sub-supercritical transition occurring over the second bump at x = x2 , the inverse transition is revealed by a sudden decrease of the Froude number, caused by an abrupt increase of the depth and a corresponding decrease of the flow velocity, i.e. a bore or hydraulic jump occurring at 31.3 ≤ x/σ ≤ 31.4. It is worth noting that the bed profile has been chosen in such a way to make this super-subcritical transition occur at a location where the bed slope, though very small, is not null: the analytical solution for both the depth and the flow velocity profile immediately downstream the hydraulic jump thus consists in a smooth variation toward the downstream horizontal profile. The simulation, though intrinsically 1D, has been performed in a 2D domain discretized by 1000 × 5 nodes. Free-slip conditions have been imposed on all boundaries. For this test case, the parameters a, b, c 0 have been set to 1, 2, gh0 respectively. A set of 81 velocities has been employed. The relaxation time was chosen as τ ∗ = 0.8. Courant number was set to 0.1. This test and the previous one demonstrate how the proposed model is able to flawlessly simulate 1D supercritical shallow water regimes with Froude number much higher than one. It is interesting to observe that the numerical steady flow agrees remarkably well with the analytical solution, which has been obtained in absence of dissipative forces. This suggests that the main source of numerical viscosity for the model under consideration stems from the first order finite difference discretization, and not from the viscosity of the underlying kinetic model. This will be shown in detail in Section 6.
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
125
Fig. 3. 1D steady flow profiles over Gaussian bumps, comparison between analytical (solid traces) and numerical (markers) results: upper panel, dimensionless water depth; lower panel, Froude number.
5.3. 2D circular dam-break The 2D circular dam break test of Billett and Toro [4] is here considered. The problem consists in the collapse of a circular cylinder of liquid, of initial depth H and radius R, surrounded by fluid at rest with initial depth h. The radial symmetry can be exploited in order to numerically compute a 1D shallow water reference numerical solution, which depends only on the radial coordinate and time: here a classical Runge Kutta 22 method has been used to integrate the 1D set of equations; the reference solution for the case h/ H = 1/10, R / H = 3.5 was generated discretizing R with 200 computational points, and used to assess DBM accuracy. In Fig. 4 the water depth and Froude number are shown at two different instants of time. Froude number reaches Fr 2 during the simulation, making this test ideal for assessing the capability to simulate 2D transcritical flows. The considered simulation has been run on a 600 × 600 Cartesian grid, Courant number was set to 0.1. From Fig. 4 it is evident that the DBM numerical solution, which has been formulated as 2D, matches remarkably well the reference one. The smearing visible at some locations is imputable to the numerical viscosity and will be discussed in Section 6. Error has been calculated for the Froude number and water depth using a formula similar to (29) but accounting
126
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
Fig. 4. Circular dam break: comparison in terms of Froude number (dashed lines and empty markers, left panel), and water depth (continuous lines and filled markers, right panel), for two instants (gray traces t = 0.2, black traces t = 0.5).
for problem axial symmetry. A value slightly lower than 10−2 for the depth and 10−1 for the Froude number has been found for the case considered. 5.4. 2D dam-break The 2D dam break test of Fennema and Chaudhry [6] is here considered. It consists of the propagation of a wave triggered by the instantaneous collapse of a lock separating two parts of the domain, each one with a different initial water level at rest (see the sketch in Fig. 5 for reference). The square domain has a side of 200 m, is divided into two equal parts by a 10 m thick wall. The wall has a 75 m wide breach, extending from y = 95 m to y = 170 m. The initial water level is equal to hm on the left of the wall (namely x ≤ 95 m), and to h v on the right of the wall (namely x > 95 m). At t = 0 s the breach is considered open and a 2D dam break flow is generated. Numerical simulations described in the following are relative to a case with hm = 10 m, h v = 0.3 m. Such markedly high difference was chosen to produce transcritical 2D shocks. The simulation is characterized by the propagation of a weak, rounded-shaped shock traveling almost perpendicularly to the dividing wall, meanwhile, a rarefaction wave spreads radially into the left part of the domain with a speed which is almost half of the shock’s one ((27) gives a close estimate of such celerities). The shock impacts on the opposite wall and is reflected in the form of a strong, wide and backward traveling hydraulic jump; the accumulation of water along the opposite wall spreads laterally and impinges into both lateral walls creating a complicated system of waves traveling back toward the breach and propagating upwind through the aperture into the still basin. No analytical solution is available for this benchmark thus in this work a well established and validated numerical model has been considered as reference. The reference numerical model integrates the 2D shallow water equations (1) by means of a Finite Volume, shock capturing scheme over an unstructured triangular mesh [14]. It employs a second order Total Variation Diminishing, Weighted Average Flux method [23]. Such scheme has been developed in order to guarantee correct propagation speed of discontinuous solutions, while maintaining a second order accuracy over smooth ones. The same benchmark has been simulated also by means of the GK scheme [16] in order to compare the accuracy and the efficiency of the proposed multispeed approach with another Boltzmann-based algorithm able to simulate supercritical flows. The DBM and the GK simulations have been carried out on a 100 × 100 grid. Values of the DBM parameters a, b, c 0 have been set to 1, 2, ghm respectively. The employed velocity set retained 81 elements and τ ∗ was set to 0.8. Courant number has been set to 0.1. The reference numerical simulation has been carried out on an unstructured mesh of 3574 triangular cells. The assessment of the accuracy of the present model is carried out comparing the time histories of depth and Froude number (Fig. 6), and specific discharges q x , q y (Fig. 7) at points P 1 , P 2 , P 3 , P 4 (see Fig. 5 for points’ locations), with the results yielded by the reference model. Fig. 6 and Fig. 7 include GK results. Water depths are scaled with hm . Specific discharges q x , q y are defined as: q x = uh, q y = vh and scaled with q0 = hm
ghm . Finally, time is scaled with t 0 = hm . ghm
The three models substantially yield the same results, being the shape of the time histories very similar at most of the measuring points. The propagation of steep shocks, observable at P 2 for example, is well reproduced both in terms of strength and speed, with a slight advance of the proposed model for what concerns the bore reflected by the front wall. It is worth noting that the flow becomes markedly supercritical, as can be seen by inspecting the Froude number values attained at points P 1 and P 2 in Fig. 6. The proposed model shows a marked numerical viscosity compared to the reference
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
127
Fig. 5. Sketch of the 2D dam break benchmark of Fennema and Chaudhry [6].
Fig. 6. Time histories of dimensionless water depth (black trace) and Froude number (gray trace) at gauge locations. Comparison among DBM numerical model (solid traces), GK numerical model (dotted traces) and reference solution (dashed traces).
and the GK model, due to both the adopted time–space integration of the governing (26) and the value of the relaxation time. Table 1 shows the values of mean absolute errors Err [h], Err[Fr] (29) between the two models (the DBM and the GK) and the reference one in order to compare the numerical accuracy. The higher accuracy of the proposed DBM with respect to that of the GK model is evident. The main purpose of this benchmark is, as stated above, to check for the ability of the model to correctly compose multidirectional transcritical shocks: based on the outcome of the comparison, it can be concluded that the model possesses such feature. In addition to such quantitative assessment, a general idea on the ability
128
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
Fig. 7. Time histories of dimensionless specific discharge (along x-direction black traces, along y-direction, gray traces) at gauge locations. Comparison among DBM numerical model (solid traces), GK numerical model (dotted traces) and reference solution (dashed traces).
Table 1 Comparison between the mean absolute error attained by DBM and GK numerical results with respect to the reference solution at points P 1 , P 2 , P 3 and P 4 .
P1 P2 P3 P4
DBM–REF Err (h)
GK–REF Err (h)
DBM–REF Err(Fr)
GK–REF Err(Fr)
0.06 0.06 0.03 0.02
0.09 0.06 0.05 0.04
0.11 0.09 0.12 0.26
0.31 0.12 0.26 0.24
Fig. 8. 2D plot of the water depth in the whole domain at t = 26.9 s. Left panel: DBM numerical results. Right panel: numerical reference results. Distances and water depth values are expressed in meters.
of the proposed multispeed DBM can be gleaned from Fig. 8, which shows the distribution of the water surface elevation in the whole domain at t = 26.9 s. In Fig. 8 both the horizontal and the vertical lengthscales are expressed in meters. All the most important flow structures are very similar in both panels, in particular the shape and the position of the curved shock. Finally, as far as the computational performance is concerned, the DBM, due to its intrinsically simple structure, presents a slighter advantage in terms of computational time with respect to both the GK and the reference model. However, further dedicated tests are needed to draw any conclusion on the computational burden needed by the proposed method.
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
129
6. Kinetic and numerical viscosity The Shallow Water Equations are usually expressed as an inviscid set of conservation laws. All the benchmarks employed in this work compare the results of the DBM with inviscid solutions (either analytical or numerical) of the system. The DBM, while showing an overall considerable agreement with those reference solutions, inevitably introduces additional viscosity. The two sources of dissipation are the viscosity term νkin introduced by the collision step in the kinetic equation (26) and the numerical diffusion νFD generated by the Finite Difference scheme employed to solve it. The latter can be expressed as in Fletcher [7]: νFD = xc 0 (1 − C )/2 where C is the Courant number imposed. The ratio between the two viscosities can be then expressed as:
νFD /νkin =
1−C
(34)
2C τ ∗
The expression for the kinetic viscosity is deduced from (21) estimating τ ∼ τ ∗ t. All simulations were run with C = 0.1 and τ ∗ = 0.8, thus resulting in the numerical viscosity being roughly five times greater than the kinetic one. 7. Conclusions In this work a Discrete Boltzmann Model able to solve the 2D transcritical shallow water equations has been developed and validated. The model employs an original discretization of the continuous particle velocity space consisting of two sets of velocities grouped on the basis of their magnitude. The particular structure of the chosen velocity set allows to significantly reduce the number of unknown coefficients of a fourth order polynomial expression of the Equilibrium Distribution Functions: the coefficients are obtained by matching discrete hydrodynamics moments up to fourth order with their continuous counterparts. Several benchmark cases have been carried out in order to allow for a thorough assessment of the ability of the proposed DBM to correctly converge to the solution of the shallow water equations when 1D or 2D supercritical flows develop. This study shows that Boltzmann-based methods can be extended to the simulation of trans- and supercritical shallow water flows, frequently found in real-world applications. Acknowledgement The research has been supported by the Italian national project “Hydroelectric energy by osmosis in coastal areas”, PRIN 2010–2011. Appendix A. Definition of the EDF coefficients
If we define φ = c 0 / gh, for k = 0:
A 0 = 1 − a2 + b 2 −
C0 =
4
φ2
/ φ 2 a2 b 2 2
2
φ
− a2 − b 2 / a2 b 2
G 0 = 1/ a 2 b 2
Bk = β
2(2 − φ 2 b2 ) a2 φ 4 4(2 − φ 2 b2 ) a2 φ 2
C k = −β
Dk = β
Ek = β
2(2 − φ 2 b2 ) a2 φ 2
8(3 − φ 2 b2 ) a4 φ 2 16 3a4
(A.2)
Defining β = 1/[ N T (a2 − b2 )], for 1 ≤ k ≤
Ak = β
(A.1)
(A.3) NT 2
:
(A.4)
(A.5)
(A.6)
(A.7)
(A.8)
130
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
Fk = 0
(A.9)
G k = −β
1
(A.10)
a2
Hk = 0 Ik = β For
NT 2
(A.11) 8
(A.12)
a6
+ 1 ≤ k ≤ NT : 2(a2 φ 2 − 2)
Ak = β
4(a2 φ 2 − 2)
Bk = β
2(a2 φ 2 − 2) b2 φ 2
8(a2 φ 2 − 3)
Dk = β
16(2a2 − 3b2 ) 3b6 8(a2 − b2 )
b4
H k = −β
24(a2 − b2 )
8 b8
(A.17)
(A.18)
b4
1 2 3a − 2b2
Gk = β
(A.15)
(A.16)
b4 φ 2
F k = −β
Ik = β
(A.14)
b2 φ 2
C k = −β
Ek = β
(A.13)
b2 φ 4
b6 3a2 − 4b2
(A.19)
(A.20)
(A.21)
Appendix B. Chapman–Enskog expansion of the kinetic equations (26) Consider the scaling:
t=
c0t L
= L∇ ∇ τ=
= ck =
τ
tc c0tc L
ck c0
fk fk = , H
fe = k
f ke H
(B.1)
where H , L, t c are respectively the macroscopic vertical and horizontal lengthscales and the mesoscopic timescale, i.e. the time interval between two successive collisions. Applying the scaling (B.1) to time and space derivatives in (26) and expanding them in perturbative series with respect to the smallness parameter (namely the Knudsen number, provided a convective scaling is assumed), the following expressions can be obtained:
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
131
fk = f e + f 1 + 2 f 2 + 3 f 3 + 4 f 4 + . . . k
k
k
k
k
∂ ∂ ∂ ∂ ∂ = + + 2 + 3 + ... ∂ t ∂ t ∂ t ∂ t4 ∂t 1 2 3 = ∇1 + ∇2 + 2 ∇3 + 3 ∇4 + . . . ∇
(B.2)
The following equations are obtained at each order:
e ∂ f ke f1 + ∇1 · ck fk = − k ∂ t1 τ e e ∂ ∂ f k1 fk f2 =− k + ∇1 · ck f k1 + + ∇2 · ck fk ∂ t1 ∂ t2 τ 1 e 2 e ∂ fk f3 ∂ fk ∂ fk =− k + ∇1 · ck f k2 + + ∇2 · ck f k1 + + ∇3 · ck fk ∂ t1 ∂ t2 ∂ t3 τ 2 1 e 3 e ∂ fk f4 ∂ fk ∂ fk ∂ fk 3 2 1 + + + ∇1 · ck f k + + ∇2 · ck f k + ∇3 · ck f k + ∇4 · ck f k = − k ∂ t1 ∂ t2 ∂ t3 ∂ t4 τ
0) 1) 2) 3) .. .
(B.3)
Summing (B.3) with respect to k and accounting for the fact that: NT
h fe = k
k =0
H
NT hu ck f ke =
c0 H
k =0 NT
f ki = 0;
i = 1, 2, . . .
k =0 NT ck f ki = 0;
i = 1, 2, . . .
(B.4)
k =0
it is straightforward to obtain the shallow water mass balance equation, which in dimensional form is given by:
∂h + ∇(hu) = 0 ∂t
(B.5)
Multiplying each (B.3) by ck , summing with respect to k and accounting for the definitions of the hydrodynamic moments:
2 =
NT h2 1 ck ck f ke = 2 hu ⊗ u + g I k =0
3 =
c0 H
2
NT ck ck ck f ke k =0
NT 4 = ck ck ck ck f ke k =0
where symbol ⊗ indicates the dyadic product of the vector u with itself and I is the identity matrix, the following form of the macroscopic momentum balance equation is obtained:
c 02 H h2 ∂ hu ∂ 2 + ∇ · hu ⊗ u + g I = τ ∇1 · + ∇1 · 3 ∂t 2 L ∂ t1 ∂ 2 ∂ 2 2 2 + ∇1 · 3 + τ ∇1 · + ∇2 · 3 + τ ∇2 · ∂ t1 ∂ t2 2 ∂ ∂ 2 3 + τ ∇1 · + 2 ∇ · + ∇ · ∇ · . . . − 2 1 3 1 1 4 ∂ t1 ∂ t 12
(B.6)
132
M. La Rocca et al. / Journal of Computational Physics 284 (2015) 117–132
Thus (26), together with the velocity set (18) and the EDFs (16) are equivalent to the shallow water momentum balance equation (2). The equivalence is meant in the limit of small . The approximation error is proportional to 3 , due to the fact that the proposed EDFs (16) together with the velocity set (18) are able to reproduce exactly all the hydrodynamic moments appearing at right hand side of (B.6): i.e. the first term which is not correctly reproduced in (B.6) is proportional to 3 . References [1] Milton Abramowitz, Irene A. Stegun, Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, Courier Dover Publications, 2012. [2] Cyrus K. Aidun, Jonathan R. Clausen, Lattice-Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (1) (2010) 439–472. [3] R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation: theory and applications, Phys. Rep. 222 (3) (1992) 145–197. [4] S.J. Billett, E.F. Toro, On WAF-type schemes for multidimensional hyperbolic conservation laws, J. Comput. Phys. 130 (1) (1997) 1–24. [5] Bastien Chopard, V.T. Pham, Laurent Lefevre, Asymmetric lattice Boltzmann model for shallow water flows, Comput. Fluids 88 (2013) 225–231. [6] Robert J. Fennema, M. Hanif Chaudhry, Explicit methods for 2-D transient free surface flows, J. Hydraul. Eng. 116 (8) (1990) 1013–1034. [7] C.A.J. Fletcher, Computational Techniques for Fluid Dynamics, Springer, Berlin, 1991. [8] M.S. Ghidaoui, J.Q. Deng, W.G. Gray, K. Xu, A Boltzmann based model for open channel flows, Int. J. Numer. Methods Fluids 35 (4) (2001) 449–494. [9] M.S. Ghidaoui, A.A. Kolyshkin, J.H. Liang, F.C. Chan, Q. Li, K. Xu, Linear and nonlinear analysis of shallow wakes, J. Fluid Mech. 548 (2006) 309–340. [10] Charles Hirsch, Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics, vol. 1, Butterworth– Heinemann, 2007. [11] Taehun Lee, Ching-Long Lin, A characteristic Galerkin method for discrete Boltzmann equation, J. Comput. Phys. 171 (1) (2001) 336–356. [12] J.H. Liang, M.S. Ghidaoui, J.Q. Deng, W.G. Gray, A Boltzmann-based finite volume algorithm for surface water flows on cells of arbitrary shapes, J. Hydraul. Res. 45 (2) (2007) 147–164. [13] Chiang C. Mei, The Applied Dynamics of Ocean Surface Waves, vol. 1, World Scientific, 1989. [14] Pietro Prestininzi, Numerical modelling of fluvial floods, PhD thesis, Università degli Studi Roma TRE, 2009. [15] Pietro Prestininzi, Michele La Rocca, Reinhard Hinkelmann, et al., Comparative study of a Boltzmann-based finite volume and a lattice Boltzmann model for shallow water flows in complex domains, Int. J. Offshore Polar Eng. 24 (03) (2014) 161–167. [16] Pietro Prestininzi, Michele La Rocca, Andrea Montessori, Giampiero Sciortino, A gas-kinetic model for 2D transcritical shallow water flows propagating over dry bed, Comput. Math. Appl. 68 (4) (2014) 439–453. [17] Pietro Prestininzi, Giampiero Sciortino, Michele La Rocca, On the effect of the intrinsic viscosity in a two-layer shallow water lattice Boltzmann model of axisymmetric density currents, J. Hydraul. Res. 51 (6) (2013) 668–680. [18] X.W. Shan, X.Y. He, Discretization of the velocity space in the solution of the Boltzmann equation, Phys. Rev. Lett. 80 (1) (1998) 65–68. [19] Xiaowen Shan, Xue-Feng Yuan, Hudong Chen, Kinetic theory representation of hydrodynamics: a way beyond the Navier–Stokes equation, J. Fluid Mech. 550 (2006) 413–441. [20] James D. Sterling, Shiyi Chen, Stability analysis of lattice Boltzmann methods, J. Comput. Phys. 123 (1) (1996) 196–206. [21] James Johnston Stoker, Water Waves: The Mathematical Theory with Applications, vol. 36, John Wiley & Sons, 2011. [22] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford Science Publications, Oxford, 2001. [23] E.F. Toro, A weighted average flux method for hyperbolic conservation laws, Proc. R. Soc. Lond. A 423 (1865) (1989) 401–418. [24] S. Ubertini, G. Bella, S. Succi, Lattice Boltzmann method on unstructured grids: further developments, Phys. Rev. E 68 (1) (2003) 1–10. [25] Stefano Ubertini, G. Bella, Sauro Succi, Unstructured lattice Boltzmann equation with memory, Math. Comput. Simul. 72 (2) (2006) 237–241. [26] Pham Van Thang, Bastien Chopard, Laurent Lefèvre, Diemer Anda Ondo, Eduardo Mendes, Study of the 1D lattice Boltzmann shallow water equation and its coupling to build a canal network, J. Comput. Phys. 229 (19) (2010) 7373–7400. [27] Minoru Watari, Michihisa Tsutahara, Possibility of constructing a multispeed Bhatnagar–Gross–Krook thermal model of the lattice Boltzmann method, Phys. Rev. E 70 (1) (2004) 016703. [28] Kun Xu, A well-balanced gas-kinetic scheme for the shallow-water equations with source terms, J. Comput. Phys. 178 (2) (2002) 533–562. [29] Jian Zhou, Lattice Boltzmann Methods for Shallow Water Flows, Springer, Berlin, 2004.
