Application of a Steady Meandering River with Piers Using a Lattice Boltzmann Sub-Grid Model in Curvilinear Coordinate Grid

: A sub-grid multiple relaxation time (MRT) lattice Boltzmann model with curvilinear coordinates is applied to simulate an artiﬁcial meandering river. The method is based on the D2Q9 model and standard Smagorinsky sub-grid scale (SGS) model is introduced to simulate meandering ﬂows. The interpolation supplemented lattice Boltzmann method (ISLBM) and the non-equilibrium extrapolation method are used for second-order accuracy and boundary conditions. The proposed model was validated by a meandering channel with a 180 ◦ bend and applied to a steady curved river with piers. Excellent agreement between the simulated results and previous computational and experimental data was found, showing that MRT-LBM (MRT lattice Boltzmann method) coupled with a Smagorinsky sub-grid scale (SGS) model in a curvilinear coordinates grid is capable of simulating practical meandering ﬂows.


Introduction
Derived from the Lattice Gas Automata (LGA) [1,2], the single relaxation method (called the lattice Bhatnagar Gross Krook (LBGK) method) [3,4] is a promising and powerful tool for computational fluid dynamics. It has also been successfully applied to simulate various of flow problems, such as free surface flow [5], advection and dispersion problems [6], multiphase fluids [7,8], and shallow water flows [9,10]. However, numerical instability is one of problems of the LBGK method, especially when the Reynolds number is high or the viscosity is low [11,12]. The multiple relaxation time (MRT) lattice Boltzmann method was proposed and developed [13][14][15] to overcome these shortcomings; by establishing a model on moment space rather than on discrete space, different relaxation times can be chosen for different moments, which leads to an improvement in the stability of the LBGK method [15].
Standard lattice Boltzmann method (LBM) is restricted to regular lattices, which is limits the simulation of curved and complex natural rivers. One way to solve this problem is to use nonuniform lattices. In recent years, different methods have been developed to extend the LBM on a nonuniform mesh, including the interpolation-supplemented scheme (ISLBE) [16], grid refinement scheme [17][18][19], dynamically adaptive grids for shallow water simulations [20], and the MRT-LBM for transformed equations in a curvilinear coordinates system [21].
The interpolation-supplemented scheme (ISLBE) was first proposed to employ nonuniform rectangular grids [16], and it was successfully applied to the flow around a circular cylinder in a polar coordinate grid system under different Reynolds numbers [22]. Shyam Sunder et al. [23] investigated the parallel performance of the ISLBE scheme and demonstrated that the ISLBE scheme could obtain a good parallel performance, although it increased the communication and computational time, Moreover, the generalized form of the interpolation supplemented lattice Boltzmann method (GILBM) was proposed to simulate the steady flow in generalized coordinate [24]. Qu et al. [25] applied the isoparametric transformation to the ISLBE and therefore arbitrarily structural grids could be used. More recently, Zhao applied the GILBM to shallow water equations, allowing the flow problem in curved and meandering open channels to be accurately resolved based on a curvilinear coordinate grid system [26].
The phenomenon of bridge piers in a river is a common flow problem. On account of the problem of piers holding up water and the phenomenon of flow around cylinders, a suitable turbulence model is significant. In recent years, large eddy simulation (LES) has been widely used and studied for turbulent flow. Hou et al. [27] first integrated a sub-grid scale (SGS) stress model with LBGK for turbulence modeling. Zhou [28] extended the scheme to shallow water flows. Yu et al. [29] integrated MRT-LBM with LES and demonstrated the method's superiority to LBGK-LES. However, curved channels are seldom considered in these studies. Therefore, in this study we attempt to introduce a standard Smagorinsky sub-grid scale (SGS) model into Zhao's scheme to simulate curved steady flows with piers. In addition, a 180 • open channel is used to validate the model.
The rest of this study is organized as follows. The governing equations of shallow water flow and a sub-grid lattice Boltzmann method with curvilinear coordinates are described in Section 2. Section 3 presents simulating results and discussion. Conclusions are summarized in Section 4.

Governing Equations
Shallow water equations were used for our simulation, which can be written as: In Equations (1) and (2), h is water depth; x i is the distance in the i direction; u j is depth-averaged velocity components in the j direction; t is the time, g = 9.81 m/s 2 is the gravitational acceleration; ν t = ν k + ν e is the total viscosity.
The kinematic viscosity ν k is usually defined as in which e = δx/δt, δx is the lattice size; δt is the time step. In the Smagorinsky model [30], the eddy viscosity ν e is given by where C s is the Smagorinsky constant (C s = 0.15 in the present study), l s is the length scale, and S ij is the strain rate tensor given by [27,28] S ij = 1 2h F i is the force term; without the Coriolis term, it can be defined as: where ρ is the water density, z b is the bed elevation, and τ bi is the bed shear stress: Here, C b is the bed friction coefficient and C b = g/C 2 z , where C z is the Chezy constant given by the Manning equation, C z = h 1/6 /n b , and n b is the Manning coefficient.

A Sub-Grid Lattice Boltzmann Model
It is simple to introduce SGS into MRT-LBM. By adding the calculation module of eddy relaxation time to the whole calculation, a sub-grid lattice Boltzmann model can be established.
In this study, MRT-LBM is used to solve shallow water equations, and it has already been applied successfully by numerous researchers [6,29,31,32]. Our simulation is based on the D 2 Q 9 model and space is discretized into a nine-speed square lattice (see Figure 1). The particle velocity e α is defined as: Water 2018, 10, x FOR PEER REVIEW  3 of 14 where is the water density, is the bed elevation, and is the bed shear stress: Here, is the bed friction coefficient and = / , where is the Chezy constant given by the Manning equation, = ℎ ⁄ ⁄ , and is the Manning coefficient.

A Sub-Grid Lattice Boltzmann Model
It is simple to introduce SGS into MRT-LBM. By adding the calculation module of eddy relaxation time to the whole calculation, a sub-grid lattice Boltzmann model can be established.
In this study, MRT-LBM is used to solve shallow water equations, and it has already been applied successfully by numerous researchers [6,29,31,32]. Our simulation is based on the model and space is discretized into a nine-speed square lattice (see Figure 1). The particle velocity is defined as: where is the particle distribution function; the relationship between the distribution function and the moment is = , = and the bold-face symbols denote nine-dimensional column vectors, e.g., = , , ⋯ ; is the equilibrium distribution functions of the moment ; is the post-collision state. Ginzburg [13] first proposed the general form of the transformation matrix which can be defined as [6,15,31]: is associated with the force term. In a model, the corresponding equilibrium distribution functions in the moment space are expressed as [6,15,31]: The MRT-LBM contains two steps: a collision step and a streaming step:

Collision and Forcing
where f α is the particle distribution function; the relationship between the distribution function and the moment isf = Mf, f = M −1f and the bold-face symbols denote nine-dimensional column vectors, e.g., f = [ f 0 , f 1 , · · · f 8 ] T ;f eq α is the equilibrium distribution functions of the momentf α ; f α is the post-collision state. Ginzburg [13] first proposed the general form of the transformation matrix M which can be defined as [6,15,31]: F i is associated with the force term. In a D 2 Q 9 model, the corresponding equilibrium distribution functions in the moment spacef eq α are expressed as [6,15,31]: S is the relaxation matrix in the moment space: where With the eddy viscosity ν e of the SGS model in consideration, s 8 and s 9 are decided as where τ is the single relaxation time, τ t is the total relaxation, and τ e is the eddy relaxation time, which is given by [27,28]: in which e αi e αj f α − f eq α . (16) e αi is the velocity vector of a particle in the i spatial coordinate. In LBM, the equilibrium distribution functions (EDFs) can be calculated by a Taylor series expansion of a Maxell-Boltzmann equilibrium distribution. According to constraint conditions of EDFs, the local equilibrium distribution f eq α for shallow water equations can be computed by the method of undefined coefficients [9,33], and f eq α can be expressed as: where ω α is the weight coefficient, ω 0 = 4/9, ω 1−4 = 1/9, ω 5−8 = 1/36. The water depth h and velocity u can be calculated by the distribution function below:

Curvilinear Coordinates
The lattice Boltzmann model on curvilinear coordinates was presented by using the GILBM [24].

Collision and Forcing
in which where e α ≡ e α,i is the particle velocity of curvilinear coordinates (ε, η) and can be calculated as To calculate the contravariant velocities at each node, the estimation of the metrics is given by in which J is described as To integrate the particle velocity, the second-order two-step Runge-Kutta method is used: Second step : The interpolation function is necessary to calculate the right-hand side of Equations (20) and (27). In both equations, the value between the grid points is required.
The second-order upwind quadratic interpolation is used, and can be expressed as where g α,l,m is the stencil depending on the position of ε − ∆ε α . The coefficients are described as

Boundary Conditions
Boundary conditions are significant and can affect the accuracy of the lattice Boltzmann method. The non-equilibrium extrapolation method [34] of second-order accuracy was chosen to determine the distribution functions at the boundaries from the given macroscopic variables: where the subscript w denotes the wall nodes and f represents fluid nodes. The fluid water depth h f and macroscopic velocities u f can be computed from discharge or water level according to real cases.
For the wall nodes, since the water depth h w and the macroscopic velocities u w are not explicit, they are estimated from the neighboring fluid nodes. For the macroscopic velocities u w , they can be estimated from the neighboring fluid velocities at the slip wall and are equal to zero at the non-slip wall. For the water depth, a three-point Lagrangian formula is applied [35]: In our simulation, the criterion of steady states is defined as

Model Simulation and Discussion
In this part, an open channel flow with a 180 • bend is investigated to test the accuracy of the proposed scheme. Moreover, a real meandering river with piers is simulated to test the application of the coupled model. The results of the simulation and a discussion are presented as follows.

Open Channel Flow with a 180 • Bend
De Vriend [36] experimentally studied the 180 • open channel, and it has also been used by many researchers to validate their models [26,37,38]. Zhang [38] has established a 3D Re-Normalisation Group (RNG) k − ε turbulence model with curvilinear coordinates to simulate meandering rivers and channels. Therefore, it is appropriate to compare these results with ours. The studied channel is 1.7 m wide and the centerline radius is R = 4.25 m. There are two 6-m long straight reaches connected to the bend. The channel boundaries are hydraulically smooth and the bed slop is zero. The upstream flow discharge is 0.19 m 3 /s and the downstream value is given with the terminal water depth H 0 = 0.18 m. A uniform mesh of 100 × 26 ( Figure 2) with δx = δy = 0.047 m is employed, the particle velocity e is 2 m/s, and the Froude number Fr = 0.47. Also, the chosen time step is δt = δx/e = 0.0235 s and the relaxation time is τ = 0.819 s, as determined by Equation (3).
In our simulation, the criterion of steady states is defined as

Model Simulation and Discussion
In this part, an open channel flow with a 180° bend is investigated to test the accuracy of the proposed scheme. Moreover, a real meandering river with piers is simulated to test the application of the coupled model. The results of the simulation and a discussion are presented as follows.

Open Channel Flow with a 180° Bend
De Vriend [36] experimentally studied the 180° open channel, and it has also been used by many researchers to validate their models [26,37,38]. Zhang [38] has established a 3D Re-Normalisation Group (RNG) k − ε turbulence model with curvilinear coordinates to simulate meandering rivers and channels. Therefore, it is appropriate to compare these results with ours. The studied channel is 1.  Comparison of the water depth at the central line as well as the inner and outer banks are depicted in Figure 3. The error analysis of different models is presented in Table 1. For LBM, the maximum root mean squared error(RMSE) and the relative RMS error (RRE), which occurred at the outer bank, are 0.013% and 16.7%, respectively. In Zhang's 3D model, these values are 0.012% and 16.5%, which occurred at the center line. Simulations at the inner bank for the two models are better than that at the central line and the outer bank, since the relative RMS error at the inner bank is the smallest. In general, compared to the experimental data, our 2D scheme achieved an acceptable result in the water depth, but a closer agreement between the experimental data and results of the 3D model was realized. Comparison of the water depth at the central line as well as the inner and outer banks are depicted in Figure 3. The error analysis of different models is presented in Table 1. For LBM, the maximum root mean squared error(RMSE) and the relative RMS error (RRE), which occurred at the outer bank, are 0.013% and 16.7%, respectively. In Zhang's 3D model, these values are 0.012% and 16.5%, which occurred at the center line. Simulations at the inner bank for the two models are better than that at the central line and the outer bank, since the relative RMS error at the inner bank is the smallest.
In general, compared to the experimental data, our 2D scheme achieved an acceptable result in the water depth, but a closer agreement between the experimental data and results of the 3D model was realized.  Comparisons between predicted and measured depth-averaged velocities at six cross-sections are depicted in Figure 4; the error analyses of different sections are shown in Table 2. At section 60° and 180°, our scheme achieved better results and the RMS error values were 0.06 and 0.017, while in the remaining sections the 3D model performed better. For both methods, the maximum RRS error occurred at section 90° and the values were 0.063 and 0.058 (see Table 2). Generally, although there are some discrepancies at section 90°, the comparison of results between the measured data and predicted data was acceptable.
The results of water depth and velocities indicate that Zhang's model achieved better results. The main reason for this is that Zhang's model is established as a 3D RNG k-ε turbulence model, which considers the vertical direction and the influence of the secondary flow in the meandering channel. Nevertheless, as our scheme is a 2D model, it is simpler in programming and the results show reasonable accuracy.   (1) MAE is mean absolute error; RMSE is root mean squared error, which is the average of the squared differences between measured and predicted values; RRE is defined as the ratio of the RMS error to the observed change.
Comparisons between predicted and measured depth-averaged velocities at six cross-sections are depicted in Figure 4; the error analyses of different sections are shown in Table 2. At section 60 • and 180 • , our scheme achieved better results and the RMS error values were 0.06 and 0.017, while in the remaining sections the 3D model performed better. For both methods, the maximum RRS error occurred at section 90 • and the values were 0.063 and 0.058 (see Table 2). Generally, although there are some discrepancies at section 90 • , the comparison of results between the measured data and predicted data was acceptable.
The results of water depth and velocities indicate that Zhang's model achieved better results. The main reason for this is that Zhang's model is established as a 3D RNG k-ε turbulence model, which considers the vertical direction and the influence of the secondary flow in the meandering channel. Nevertheless, as our scheme is a 2D model, it is simpler in programming and the results show reasonable accuracy.

Meandering River with Piers
In order to test the capability of simulating the practical problem in real rivers, the present scheme was used to simulate the meandering flow of the Longhua River in Shenzhen, China. This study area is 606 m long and 31.5 m wide with two bend sections (section 1# and section 2#). Figure 5 shows a schematic description of the simulation domain. White circles depicted in the figure represent the piers in this river, across which there is a small bridge.

Meandering River with Piers
In order to test the capability of simulating the practical problem in real rivers, the present scheme was used to simulate the meandering flow of the Longhua River in Shenzhen, China. This study area is 606 m long and 31.5 m wide with two bend sections (section 1# and section 2#). Figure 5 shows a schematic description of the simulation domain. White circles depicted in the figure represent the piers in this river, across which there is a small bridge. The accuracy and stability of the traditional finite volume method (FVM) has been demonstrated in the simulation of practical rivers [39,40]. Therefore, our results were compared with FVM results and monitoring data. The monitoring data were obtained by monitoring the river on 8 August, 2017, for which water depth and depth-average velocities were monitored at the bends. There were four monitoring points set every 10 meters from the inner bank at each bend.
The main parameters are described in Table 3. The upstream discharge was 494 m 3 /s and the downstream water depth was 5.22 m. The Manning's coefficient of the bed was 0.025 and constant particle velocity was e = 10 m/s. The time step was δ = 0.05s and the single relaxation time was τ = 0.56 s. In LBM simulations, body-fitted coordinate grids were used and a uniform mesh of 316 × 20 (see Figure 6) was applied. The minimum lattice length was 1.00035 m, while the diameter of piers was 1.2 m. Like the wall boundaries, these piers were simulated as obstacles, and a slip boundary transformed on curvilinear coordinates was used.
In the FVM model, the Reynolds-averaged Navier-Stokes equations (RANS) were solved to simulate piers, and triangle meshes were used. The minimum element area was 0.66 m 2 . The accuracy and stability of the traditional finite volume method (FVM) has been demonstrated in the simulation of practical rivers [39,40]. Therefore, our results were compared with FVM results and monitoring data. The monitoring data were obtained by monitoring the river on 8 August 2017, for which water depth and depth-average velocities were monitored at the bends. There were four monitoring points set every 10 m from the inner bank at each bend.
The main parameters are described in Table 3. The upstream discharge was 494 m 3 /s and the downstream water depth was 5.22 m. The Manning's coefficient of the bed n b was 0.025 and constant particle velocity was e = 10 m/s. The time step was δt = 0.05 s and the single relaxation time was τ = 0.56 s. Notes: Triangle meshes are used in FVM, therefore it shows different meshes. This river is an artificial river, so that the bed slope is constant. R e = 156.83 is the Reynolds number, which can be defined as: R e = hU 0 /2υ, where h is the height of the entry section, υ is the kinematic viscosity, and U 0 is the maximum inlet velocity.
In LBM simulations, body-fitted coordinate grids were used and a uniform mesh of 316 × 20 (see Figure 6) was applied. The minimum lattice length was 1.00035 m, while the diameter of piers was 1.2 m. Like the wall boundaries, these piers were simulated as obstacles, and a slip boundary transformed on curvilinear coordinates was used.
In the FVM model, the Reynolds-averaged Navier-Stokes equations (RANS) were solved to simulate piers, and triangle meshes were used. The minimum element area was 0.66 m 2 . A comparison was made at the bends of the river. The water depth at the bends are plotted in Figure 7, and the total velocity is shown in Figure 8. The water depth of section 1# decreases from the outer bank to the inner bank, while at section 2# the water depth decreases from the inner bank to the outer bank. The total velocity shows a decreasing trend from the inner bank to the outer bank in section 1#, while in the bend section 2# displays an increasing trend from the inner bank to the outer bank. This outcome is consistent with the actual situation; where the water depth is higher, the water flows more slowly. Table 4 shows the error analysis of the FVM and the MRT-LBM. Both methods agree well. For example, the water depth in the FVM achieves a better result, as the RMSE is 0.025 m and the relative RMS error is 12.5% at section 2#, while for MRT-LBM the RRE is 13.0% and the RMSE is 0.026 m. Generally, there are minute differences between monitored data and simulation data--this is probably due to the influence of centrifugal force which generates the secondary flow at the bend, while the 2D computation model does not consider the vertical direction. However, the accuracy is acceptable.
According to the comparison, our scheme performs well for velocity at section 1#, while at section 2# there are some discrepancies near from the inner bank and the RRE is 29.1%. In general, the FVM was superior to the proposed model. The main reason for this are: (1) the uniform mesh of the MRT-LBM is non-orthogonal, which may lead to small deviation especially at bends; (2) the governing equations should be transformed into curvilinear coordinates in the proposed model, and the finite difference approximation of the transformation matrix may lead to some discrepancies. However, our proposed model requires much fewer CPUs as well as less time for simulation. The outcome demonstrates the advantages of LBM and is consistent with relevant reports in the literature [12,15,41].  A comparison was made at the bends of the river. The water depth at the bends are plotted in Figure 7, and the total velocity is shown in Figure 8. The water depth of section 1# decreases from the outer bank to the inner bank, while at section 2# the water depth decreases from the inner bank to the outer bank. The total velocity shows a decreasing trend from the inner bank to the outer bank in section 1#, while in the bend section 2# displays an increasing trend from the inner bank to the outer bank. This outcome is consistent with the actual situation; where the water depth is higher, the water flows more slowly. Table 4 shows the error analysis of the FVM and the MRT-LBM. Both methods agree well. For example, the water depth in the FVM achieves a better result, as the RMSE is 0.025 m and the relative RMS error is 12.5% at section 2#, while for MRT-LBM the RRE is 13.0% and the RMSE is 0.026 m. Generally, there are minute differences between monitored data and simulation data-this is probably due to the influence of centrifugal force which generates the secondary flow at the bend, while the 2D computation model does not consider the vertical direction. However, the accuracy is acceptable.
According to the comparison, our scheme performs well for velocity at section 1#, while at section 2# there are some discrepancies near from the inner bank and the RRE is 29.1%. In general, the FVM was superior to the proposed model. The main reason for this are: (1) the uniform mesh of the MRT-LBM is non-orthogonal, which may lead to small deviation especially at bends; (2) the governing equations should be transformed into curvilinear coordinates in the proposed model, and the finite difference approximation of the transformation matrix may lead to some discrepancies. However, our proposed model requires much fewer CPUs as well as less time for simulation. The outcome demonstrates the advantages of LBM and is consistent with relevant reports in the literature [12,15,41]. A comparison was made at the bends of the river. The water depth at the bends are plotted in Figure 7, and the total velocity is shown in Figure 8. The water depth of section 1# decreases from the outer bank to the inner bank, while at section 2# the water depth decreases from the inner bank to the outer bank. The total velocity shows a decreasing trend from the inner bank to the outer bank in section 1#, while in the bend section 2# displays an increasing trend from the inner bank to the outer bank. This outcome is consistent with the actual situation; where the water depth is higher, the water flows more slowly. Table 4 shows the error analysis of the FVM and the MRT-LBM. Both methods agree well. For example, the water depth in the FVM achieves a better result, as the RMSE is 0.025 m and the relative RMS error is 12.5% at section 2#, while for MRT-LBM the RRE is 13.0% and the RMSE is 0.026 m. Generally, there are minute differences between monitored data and simulation data--this is probably due to the influence of centrifugal force which generates the secondary flow at the bend, while the 2D computation model does not consider the vertical direction. However, the accuracy is acceptable.
According to the comparison, our scheme performs well for velocity at section 1#, while at section 2# there are some discrepancies near from the inner bank and the RRE is 29.1%. In general, the FVM was superior to the proposed model. The main reason for this are: (1) the uniform mesh of the MRT-LBM is non-orthogonal, which may lead to small deviation especially at bends; (2) the governing equations should be transformed into curvilinear coordinates in the proposed model, and the finite difference approximation of the transformation matrix may lead to some discrepancies. However, our proposed model requires much fewer CPUs as well as less time for simulation. The outcome demonstrates the advantages of LBM and is consistent with relevant reports in the literature [12,15,41].   Two uniform meshes of 316 × 20 and 632 × 40 (see Figure 6) were applied to investigate the flow field distribution around the piers. The results in Figure 9 show that the finer grids produce a more clear flow field around the cylinders. Figure 10 depicts a 3D visualization of the streamline and water depth, which shows small vortexes around the piers. Because of the low Reynolds number, the water flows past the piers and persists without separation.   Two uniform meshes of 316 × 20 and 632 × 40 (see Figure 6) were applied to investigate the flow field distribution around the piers. The results in Figure 9 show that the finer grids produce a more clear flow field around the cylinders. Figure 10 depicts a 3D visualization of the streamline and water depth, which shows small vortexes around the piers. Because of the low Reynolds number, the water flows past the piers and persists without separation.  Two uniform meshes of 316 × 20 and 632 × 40 (see Figure 6) were applied to investigate the flow field distribution around the piers. The results in Figure 9 show that the finer grids produce a more clear flow field around the cylinders. Figure 10 depicts a 3D visualization of the streamline and water depth, which shows small vortexes around the piers. Because of the low Reynolds number, the water flows past the piers and persists without separation.

Conclusions
In this study, a sub-grid multiple relaxation time (MRT) lattice Boltzmann model with curvilinear coordinates was developed. An open channel flow with a 180° bend was simulated to validate the model. Furthermore, the calculated results were compared with the experimental data and Zhang's results. Error analysis revealed that the 3D model was superior to the proposed 2D model; however, the 2D model only requires simple programming and its results are acceptable.
A real meandering river with piers was simulated to test the application of the scheme. The results were reliable and agreed well with those of the finite volume method (FVM). Flow field and 3D streamline with height are plotted in Figures 9 and 10, in which the flow around the piers can be clearly seen. It was shown that with a low Reynolds number, the proposed method has great potential to solve realistic problems in curved rivers.
In the future, turbulent models with a high Reynolds number can be modified based on the proposed method to solve more complex flow problems. Moreover, the advection and anisotropic dispersion equations can be combined to solve the water quality problems. These extensions will be studied in our future research.
Author Contributions: L.C. and P.H. conceived the idea; Z.Z. and L.C. built the model and programming; L.C. analyzed the data and wrote the paper. All authors reviewed the manuscript.