Numerical Analysis of Flow Phenomena in Discharge Object with Siphon Using Lattice Boltzmann Method and CFD

Abstract

This article presents numerical simulation of flow in the discharge object with the welded siphon and the free water level. The main numerical tool used in this study is the lattice Boltzmann method combined with the Volume-of-Fluid approach and the Smagorinski LES model. Some aspects of the numerical method are discussed, especially the formulation of the outlet boundary condition. The simulations are carried out with in-house software based on the open-source Palabos framework. Presented results are compared with the CFD simulations, based on the ANSYS CFX software applying the SST and SAS turbulence models and the free-surface flow modeling by means of the Volume-of-Fluid method. The evolution and interactions of main flow structures are analyzed using visualizations and the spectral analysis. All numerical simulations are verified by the experimental data obtained in the hydraulic laboratory with water circuit. A stationary flow regime has been visualized by means of PIV. Both the vertical planes and horizontal planes have been examined, focused mainly on the regions below and behind the siphon outlet. The results show a good agreement of calculated and measured complex flow structures, including time-averaged and instantaneous flow fields.

1. Introduction

The paper presents the results of numerical simulations of flows in discharge objects of pumping stations. The main objective of the paper is the investigation of the usability of the lattice Boltzmann method for the case of complex flows with free surface. The work is part of a larger project involving both numerical simulations using the standard finite volume method in the so-called Volume Of Fluid (hereinafter VOF) formulation and experimental investigations using the Particle Image Velocimetry (hereinafter PIV) method. The presented paper focuses mainly on numerical simulations of two-phase free-flow in an outlet object with a siphon. This is a very complicated case where a stream of water flows into the outlet pool through a siphon with an outlet diffuser, and interacts with still water in the outlet tank and with the tank bottom. The flow field is non-stationary and contains quite complex vortex structures.

Very few studies dealing with numerical simulations of this problem can be found in the literature. The paper [1] deals with the shape optimization of the output object. However, the authors of the mentioned paper consider only single-phase stationary flow of an incompressible fluid described by a system of Navier–Stokes equations combined with a standard two-equation turbulence model. The more advanced models considering non-stationary two-phase flow with free surface [2,3] are mainly based on the numerical solution of the Navier–Stokes equations in the VOF formulation.

One of the drawbacks of simulations based on the solution of the Navier–Stokes equations is the rather tedious process of mesh generation in the case of complex geometry. The computational mesh must be significantly refined near the free surface and, in the case of large eddy simulations, a sufficiently fine mesh is also needed inside the whole region. The construction of such a mesh can be very challenging. One possible improvement is, for example, the use of adaptive meshes with automatic refinement [4,5]. However, this introduces another problem related, among others, to the efficient use of parallel computers.

An alternative method for flow simulation is the so-called Lattice Boltzmann Method (hereinafter LBM) based on the kinetic theory originally proposed in [6]. The original method was proposed for the case of compressible isothermal low Mach flows, but it was later modified for other types of problems, see e.g., the book [7]. The article [8] describes the application for convection-diffusion equation in a 2D channel, ref. [9] applies the method to two-phase flows in porous media, paper [10] deals with the simulation of turbulent flows. Paper [11] proposes the so-called discrete velocity LBM.

In this paper, we use LBM in combination with the free-surface flow model proposed in [12] for the case of foam modeling. The method has been applied for free-surface water flows [13] and then validated for many different free-surface flows. Let us mention the application of the method combined with the Smagorinsky large eddy model presented in [13] for the case of falling droplet and for the dam break problem. Articles [14,15,16] show validation of the free-surface LBM for the case of dam break. Paper [17] presents validation of the method for flows in narrow channels, ref. [18] validates the enhanced version of the method for the case of wave breaking. As far as the authors know, there are no available studies in the literature dealing with the application of LBM for the simulation of flows in discharge objects with siphon.

There are other variants of LBM for multiphase flows. The oldest one proposed by Gunstensen in [19] is based on a coloring. The model uses two sets of distribution functions and the collision and streaming is done separately for each phase with an additional collision term due to phase interface. This model suffers from some limitations such us problems with stability, expensive recoloring, and spurious currents near the interface [7,20]. Another family are the so-called pseudo-potential models such as the Shan-Chen model [21] which calculate the interaction between fluid particles using intermolecular forces. The drawback of the model is its worse stability and spurious currents especially in the case of large density ratios. The other family of models proposed in [22] incorporates the effects of phase interface by a free energy. The original model has been extended in [23] for large density ratios and should be therefore usable for two-phase water-air simulations.

The free-surface LBM proposed in [12] has been chosen for current study thanks to its computational efficiency and good stability.

2. Materials and Methods

2.1. Formulation of the Problem

The article deals with the simulation of flow problems experimentally investigated in the hydraulic laboratory of the Centre of Hydraulic Research [2]. The test section consists of 5 $\mathrm{m}$ long channel with a sluice gate at the outlet. The width of the test channel is $0.304$ $\mathrm{m}$. The water flows into the channel through a welded siphon composed of a flange DN125 which turns into rectangular cross section $0.122$ $\mathrm{m}$× $0.103$ m and ends with a rectangular diffuser with an opening angle 10 ${}^{\circ }$, see Figure 1. See [3] for detailed description of the test circuit.

The steady volume flow rate of water is considered at the inlet of the siphon (the leftmost position in Figure 1). The desired water level in the channel was kept either by setting the sluice gate at the outlet (the rightmost position in the Figure 1) in the case of experimental methods, or with an appropriate boundary condition in the case of numerical simulations (see below).

2.2. Lattice Boltzmann Method

The motion of an incompressible fluid can be described in the framework of LBM by the following set of advection-reaction PDEs for the so-called discrete density distribution functions $f_i(\overrightarrow{x},t)$, see e.g., [7]

$$\frac{\partial f_i}{\partial t}+{\overrightarrow{c}}_i·\nabla f_i=-\frac{1}{{\tau }_c}\left(f_i-f_i^{\left(eq\right)}\right)+F_i.$$

(1)

Here ${\overrightarrow{c}}_i$ represents discrete (or lattice) velocities, ${\tau }_c$ is the relaxation time, $f_i^{\left(eq\right)}$ is the equilibrium distribution, and $F_i$ represents contributions due to the external force (i.e., gravity). Moments of the distribution functions $f_i$ are the macroscopic density and momentum

$$\begin{array}{ccc}\hfill \rho (\overrightarrow{x},t)& =& \sum _i{f}_i(\overrightarrow{x},t),\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \rho (\overrightarrow{x},t)\overrightarrow{u}(\overrightarrow{x},t)& =& \sum _i{\overrightarrow{c}}_i{f}_i(\overrightarrow{x},t).\hfill \end{array}$$

(3)

Integrating Equation (1) from time t to $t+\delta t$ along the line $\overrightarrow{x}+{\overrightarrow{c}}_{i}t$ assuming constant collision term during this time interval leads to a general formulation of LBM ($\tau ={\tau }_c/\delta t$)

$$f_i(\overrightarrow{x}+{\overrightarrow{c}}_i\delta t,t+\delta t)-f_i(\overrightarrow{x},t)=-\frac{1}{\tau }\left(f_i(\overrightarrow{x},t)-f_i^{\left(eq\right)}(\overrightarrow{x},t)\right)+\delta t{F}_i.$$

(4)

For the rest of the article, we use the so-called D3Q19 model [7] where D3 denotes three-dimensional nature of the model and Q19 means that the velocity space is discretized using 19 discrete velocities ${\overrightarrow{c}}_i$ (i.e., $i=0,...,18$) given (in lattice units) as columns of the matrix $\mathbf{C}$

$$\mathbf{C}=\left[{\overrightarrow{c}}_0|{\overrightarrow{c}}_1|...|{\overrightarrow{c}}_{18}\right]=\left[\begin{array}{ccccccccccccccccccc}\hfill 0& \hfill -1& \hfill 0& \hfill 0& \hfill -1& \hfill -1& \hfill -1& \hfill -1& \hfill 0& \hfill 0& \hfill 1& \hfill 0& \hfill 0& \hfill 1& \hfill 1& \hfill 1& \hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill -1& \hfill 0& \hfill -1& \hfill 1& \hfill 0& \hfill 0& \hfill -1& \hfill -1& \hfill 0& \hfill 1& \hfill 0& \hfill 1& \hfill -1& \hfill 0& \hfill 0& \hfill 1& \hfill 1\\ \hfill 0& \hfill 0& \hfill 0& \hfill -1& \hfill 0& \hfill 0& \hfill -1& \hfill 1& \hfill -1& \hfill 1& \hfill 0& \hfill 0& \hfill 1& \hfill 0& \hfill 0& \hfill 1& \hfill -1& \hfill 1& \hfill -1\end{array}\right].$$

(5)

Moreover, we use the nondimensional form of the lattice Boltzmann equation with $\delta t=1$ and lattice spacing $\Delta =1$, i.e.,

$$f_i(\overrightarrow{x}+{\overrightarrow{c}}_i,t+1)=f_i(\overrightarrow{x},t)-\frac{1}{\tau }\left(f_i(\overrightarrow{x},t)-f_i^{\left(eq\right)}(\overrightarrow{x},t)\right)+F_i.$$

(6)

Discrete equilibrium distribution functions are defined according to [24] as

$$f_i^{\left(eq\right)}=w_i\left[\rho +\frac{{\overrightarrow{c}}_i·\overrightarrow{u}}{c_s^2}+\frac{{\left({\overrightarrow{c}}_i·\overrightarrow{u}\right)}^2}{2c_s^4}-\frac{\overrightarrow{u}·\overrightarrow{u}}{2c_s^2}\right]$$

(7)

where $c_s$ is a fictitious speed of sound equal to $1/\sqrt{3}$ for D3Q19 model and corresponding weights $w_i$ are $w_0=1/3$, $w_i=1/18$ for $i=1,2,3,10,11,12$, and $w_i=1/36$ for $i=4,5,6,7,8,9,13,14,15,16,17,18$.

The body force $\overrightarrow{F}$ (i.e., gravity) is accounted by including a forcing term according to [25] where $F_i$ are given as proposed in [25] as

$$F_i=w_i\left[\frac{\overrightarrow{B}·{\overrightarrow{c}}_i}{c_s^2}+\frac{\overrightarrow{\overrightarrow{D}}:({\overrightarrow{c}}_i\overrightarrow{c_i}-c_s^2\overrightarrow{\overrightarrow{I}})}{2c_s^4}\right],$$

(8)

where

$$\begin{array}{ccc}\hfill \overrightarrow{B}& =& \left(1-\frac{1}{2\tau }\right)\overrightarrow{F},\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill \overrightarrow{\overrightarrow{D}}& =& \left(1-\frac{1}{2\tau }\right)(\overrightarrow{v}\overrightarrow{F}+\overrightarrow{F}\overrightarrow{v}),\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill \rho \overrightarrow{v}& =& \sum _i{\overrightarrow{c}}_i{f}_i+\frac{1}{2}\overrightarrow{F}.\hfill \end{array}$$

(11)

It can be shown with the Chapman–Enskog expansion that LBM with the D3Q19 model with forcing term and with the discrete equilibrium functions (7) approximates the system of hydrodynamic equations

$$\begin{array}{ccc}\hfill \frac{\partial \rho }{\partial t}+\nabla ·\overrightarrow{u}& =& 0,\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill \frac{\partial \overrightarrow{u}}{\partial t}+\nabla ·\left(\overrightarrow{u}\overrightarrow{u}\right)& =& -\nabla p+\nabla ·\left[\nu (\nabla \overrightarrow{u}+\nabla {\overrightarrow{u}}^T)\right]+\frac{1}{\rho }\overrightarrow{F},\hfill \end{array}$$

(13)

where the pressure is related to the density as $p=c_s^2\rho$ and the viscosity $\nu$ is related to the relaxation parameter $\tau$ as

$$\nu =\left(\tau -\frac{1}{2}\right)c_s^2.$$

(14)

One can see that LBM with the equilibrium functions (7) does not correspond to the Navier–Stokes equations for incompressible fluids. Nevertheless, the error due to the density derivative in the continuity equation remains small as far as $\left|\right|\overrightarrow{u}\left|\right|<<c_s$.

To model turbulent flows at high Reynolds numbers, the large eddy simulation (LES) approach based on the Smagorinsky model proposed in [26] is used. The molecular viscosity $\nu$ is replaced by the effective viscosity ${\nu }_{eff}$ which includes also the subgrid eddy viscosity ${\nu }_{sgs}$

$${\nu }_{eff}=\nu +{\nu }_{sgs}=c_s^2\left({\tau }_{eff}-\frac{1}{2}\right).$$

(15)

The effective relaxation time ${\tau }_{eff}$ then replaces $\tau$ in (6) and takes the value of

$${\tau }_{eff}=\frac{\tau +\sqrt{{\tau }^2+2C_{s\mathrm{m}ag}{\left(\rho c_s^4\right)}^{-1}\mathsf{\Pi }}}{2},$$

(16)

where $C_{s\mathrm{m}ag}>0$ is a model constant ($C_{s\mathrm{m}ag}=0.14$ in our calculations) and $\mathsf{\Pi }=\sqrt{{\overrightarrow{\overrightarrow{\pi }}}^{\left(neq\right)}:{\overrightarrow{\overrightarrow{\pi }}}^{\left(neq\right)}}$ with

$${\overrightarrow{\overrightarrow{\pi }}}^{\left(neq\right)}=\sum _i{\overrightarrow{c}}_i{\overrightarrow{c}}_i(f_i-f_i^{\left(eq\right)}).$$

(17)

2.2.1. Free-Surface LBM

The free-surface model proposed in [13] is used in this paper. The model works with three types of cells: fluid cells, interface cells, and empty cells. Fluid cells are separated from empty cells by a layer of interface cells. The model uses an additional scalar field $ϵ(\overrightarrow{x},t)$ representing the fluid volume fraction

$$ϵ(\overrightarrow{x},t)=m(\overrightarrow{x},t)/\rho (\overrightarrow{x},t),$$

(18)

where m stands for the mass of the liquid in the lattice cell at $\overrightarrow{x}$. The mass exchange between two interface cells is approximated as

$$\Delta m_i(\overrightarrow{x},t+1)=\left[f_i(\overrightarrow{x}+{\overrightarrow{c}}_i,t)-f_i(\overrightarrow{x},t)\right]\frac{ϵ(\overrightarrow{x}+{\overrightarrow{c}}_i,t)+ϵ(\overrightarrow{x},t)}{2}.$$

(19)

In the case of fluid-fluid or fluid-interface pairs, the mass exchange is simply

$$\Delta m_i(\overrightarrow{x},t+1)=f_i(\overrightarrow{x}+{\overrightarrow{c}}_i,t)-f_i(\overrightarrow{x},t).$$

(20)

Then the mass in a given cell is updated as

$$m(\overrightarrow{x},t+1)=m(\overrightarrow{x},t)+\sum _{i=1}^{18}\Delta m_i(\overrightarrow{x},t+1).$$

(21)

The distribution functions in the interface cells corresponding to incoming velocities from empty cells must be reconstructed using the equilibrium distribution $f_i^{\left(eq\right)}$ and the nonequilibrium part from opposite directions $\overline{i}$, i.e., using a bounce back scheme, see [13]. Hence ($\overline{ı}$ corresponds to the opposite direction to i)

$$f_i(\overrightarrow{x},t+1)=f_i^{\left(eq\right)}(\overrightarrow{x},t)+f_{\overline{ı}}(\overrightarrow{x},t)-f_{\overline{ı}}^{\left(eq\right)}(\overrightarrow{x},t).$$

(22)

Once the interface cell is completely emptied ($m\le 0$), its type is changed to empty. When the interface cell becomes fully wet, it is flagged as fluid. At the end of the procedure, all empty cells directly adjacent to fluid cells are flagged as interface.

2.2.2. Boundary Conditions for Free-Surface LBM

The domain is approximated with a staircase mesh which allows simple bounce back boundary conditions at the walls, see [7]. Although the staircase mesh impairs the overall accuracy of the method, the current study focuses on the flow structures inside the outlet tank with straight walls. Therefore, the bounce back seems to be appropriate in this case.

The volumetric flow rate is given at the inlet ($x=x_{\mathrm{m}in}$) and the inlet velocity is calculated using the parabolic profile. The density at the inlet is evaluated according to [27]. Populations are divided into incoming ($i^{+}=\{10,13,14,15,16\}$), tangential ($i^0=\{0,2,3,8,9,11,12,17,18\}$), and outgoing ($i^{-}=\{1,4,5,6,7\}$). The density at the inlet ${\rho }_{in}$ is then

$${\rho }_{in}=\sum _i{f}_i=\sum _{i\in i^{+}}{f}_i+\sum _{i\in i^0}{f}_i+\sum _{i\in i^{-}}{f}_i={\rho }^{+}+{\rho }^0+{\rho }^{-}.$$

(23)

Similarly, the x component of the momentum at the inlet is

$${\rho }_{in}{u}_{in}=\sum _i{\overrightarrow{c}}_i·{\overrightarrow{e}}_1{f}_i=\sum _{i\in i+}{f}_i-\sum _{i\in i^{-}}{f}_i={\rho }^{+}-{\rho }^{-}$$

(24)

where ${\overrightarrow{e}}_1=[1,0,0]$. Eliminating unknown ${\rho }^{+}$ gives

$${\rho }_{in}=\frac{2{\rho }^{-}+{\rho }^0}{1-u_{in}}.$$

(25)

The equilibrium distribution functions are then calculated using ${\rho }_{in}$ and ${\overrightarrow{u}}_{in}=u_{in}{\overrightarrow{e}}_1$ and the nonequilibrium bounce back scheme is applied for $i\in i^{+}$

$$f_i=f_i^{\left(eq\right)}({\rho }_{in},{\overrightarrow{u}}_{in})+f_{\overline{ı}}-f_{\overline{ı}}^{\left(eq\right)}({\rho }_{in},{\overrightarrow{u}}_{in}).$$

(26)

The tangential populations $i\in {\mathsf{\beta }}^0$ are set to equilibrium distributions $f_i=f_i^{\left(eq\right)}({\rho }_{in},{\overrightarrow{u}}_{in})$.

Similar procedure is taken at the outlet. The outlet velocity is calculated using the inlet flow rate $q_{in}$, outlet channel width W, and desired water depth at the outlet H as $u_{out}=q_{in}/\left(WH\right)$. Then the water density is calculated as

$${\rho }_{out}=\frac{2{\rho }^{+}+{\rho }^0}{1+u_{out}},$$

(27)

and the nonequilibrium bounce back procedure similar to (26) is applied for $i\in i^{-}$. Finally, the volume fraction $ϵ$ and the boundary cell type are set to values from the neighbor fluid cell.

2.2.3. Computational Mesh and Case Setup

The problem was solved using a Cartesian mesh with uniform spacing $\Delta x=\Delta y=\Delta z=2\mathrm{m}\mathrm{m}$ and with a constant time step $\Delta t=5\times {10}^{-5}\mathrm{s}$. Only first $2.7\mathrm{m}$ of the channel was considered and the thickness of the siphon walls increased to $5\mathrm{m}\mathrm{m}$. The resulting mesh consists of 78 × 10⁶ cells. Figure 2 shows the side view of the computational domain. Front part of the computational domain with the siphon and part of the test channel is given in Figure 3. The test channel starts at $x=0\mathrm{m}$ with the bottom at $y=0\mathrm{m}$. The domain is symmetric around the plane $z=0\mathrm{m}$.

The use of the uniform Cartesian mesh does not allow the near wall refinement to keep $y^{+}$ below one and therefore a suitable wall model or wall function should be used. We are aware of this fact, yet we did not use any such near wall treatment in this study.

The current model solves the discrete Boltzmann equation only in fluid or interface cells (i.e., cells with non-zero fluid volume fraction) and completely ignores cells fully filled with air. Therefore, without an additional bubble model such as e.g., [28] the LBM cannot reliably predict the transient phase of the siphon fill-in. It can be seen in Figure 4 where the calculation has been started from the initial condition (28)–(30), i.e., with some empty cells inside the siphon. One should note that these “bubbles” are completely non-physical. The current version of free-surface LBM does not solve any equations inside the empty cells and simply sets the air pressure to a prescribed constant.

$$\begin{array}{ccc}\hfill \rho (\overrightarrow{x},0)& =& 1,\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill \overrightarrow{u}(\overrightarrow{x},0)& =& \overrightarrow{0},\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill ϵ(\overrightarrow{x},0)& =& \left\{\begin{array}{cc}0\hfill & \mathrm{for}>0.255\mathrm{m},\hfill \\ 1\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$

(30)

To avoid the problem with these bubbles, the calculations have been initialized with a filled-up siphon as given in Equations (31)–(33).

$$\begin{array}{ccc}\hfill \rho (\overrightarrow{x},0)& =& 1,\hfill \end{array}$$

(31)

$$\begin{array}{ccc}\hfill \overrightarrow{u}(\overrightarrow{x},0)& =& \overrightarrow{0},\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill ϵ(\overrightarrow{x},0)& =& \left\{\begin{array}{cc}0\hfill & \mathrm{for}x>0.8\mathrm{m}\wedge y>0.255\mathrm{m},\hfill \\ 1\hfill & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$

(33)

In this case, the siphon remains filled with water, see Figure 5. All calculations in the following text are therefore performed with the initial condition (31)–(33).

Distribution functions $f_i$ are then initialized to its equilibrium value. Parabolic velocity profile corresponding to the given flow rate together with water volume fraction equal to 1 is prescribed at the inlet. The inlet density is then calculated using Equation (25). The outlet velocity and density are calculated from the inlet flow rate and desired water depth as in Equation (27).

All calculations with LBM were carried out with in-house software based on an open-source CFD framework Palabos [29].

2.3. Finite Volume Method

The LBM calculations are compared with the results obtained with the standard Finite Volume Method (hereinafter FVM) in VOF formulation. The FVM calculations were carried out with the ANSYS CFX software considering a nonhomogeneous multiphase model with different velocities for the water and air fractions. The unstructured mesh with approximately 10 × 10⁶ nodes with appropriate near wall refinement was used. The maximum element size was about 3 $\mathrm{m}$$\mathrm{m}$ and the mesh was progressively refined near the walls to keep the dimensionless wall normal cell size $y_1^{+}$ bellow 1. The FVM calculations were carried out using the Scale Adaptive Simulation (SAS) model [30] and URANS with the standard two-equation SST turbulence model [31]. The fully developed velocity profile with water volume fraction equal to 1 was used as the inlet condition and a given static pressure profile was prescribed at the outlet. The dimensions of the computational domain correspond to the real geometry with an additional straight pipe of length $0.75$ $\mathrm{m}$ connected to the inlet flange. For details see [2].

2.4. Experimental Methods

The experimental data were acquired with the laser particle image velocimetry (PIV) method with one-camera arrangement of the Dantec Dynamics measurement system. Measurements were performed in 5 vertical planes located at 0 mm, $\pm 75$ mm, and $\pm 141$ mm from the vertical symmetry plane and in three horizontal planes located at 10 mm, 135 mm, and 170 mm above channel floor. The instantaneous results were captured with the frequencies of 100Hz–250Hz and the time-averaged results have used frequencies in the range 10Hz–20Hz. The relative accuracy of the measurements can be estimated as 1% to 2% in the high velocity regions. The accuracy is expected to be worse in the low velocity regions. More details on the experimental device can be found in reference [2].

3. Results

This section presents selected results of LBM and finite volume-based CFD simulations and compares these results with the experimental data obtained with the PIV measurements. The attention is focused mainly on the performance of LBM. More detailed comparison of finite volume CFD simulations with PIV measurements can be found in [2]. Both the CFD simulation and the PIV measurements were performed for two flow rates q = 13.8 L s${}^{-1}$ and q = 17.2 L s${}^{-1}$. The topology of the flow fields was in both cases very similar to the difference only in the magnitude of the velocity. Therefore, the LBM simulation was carried out only for the flow rate q = 17.2 L s${}^{-1}$ with the outlet water depth $H=268\mathrm{mm}$. Figure 6 shows the water level at the time $t=20s$ and Figure 7 the magnitude of instantaneous 2D velocity $M=\sqrt{u_x^2+u_y^2}$ in the vertical plane $z=0\mathrm{mm}$.

The unsteady results obtained with LBM have been time-averaged over the last 10 s of simulation time. The time-averaged velocity field is then compared with the time-averaged results obtained with the finite volume method and with the experiment. The results are presented in the form of the magnitude of time-averaged 2D velocity (i.e., $M=\sqrt{u_x^2+u_y^2}$ for vertical cuts and $M=\sqrt{u_x^2+u_z^2}$ for horizontal cuts) and as integral curves of the 2D velocity vectors.

Figure 8 and Figure 9 show the flow field obtained with the lattice Boltzmann method and with the standard finite volume method in VOF formulation using the SAS model. Comparison with the PIV measurements depicted in Figure 10 shows that LBM captures very well the position of the high velocity jet near the upper side of the siphon. The finite volume method overpredicts slightly the width of the jet in that case. One can clearly see that both LBM and FVM capture the dominant vortex near the bottom at $x\approx 0.35\mathrm{m}$ which can be also seen in the PIV results. Both numerical methods also predict a vortex structure above the jet ($x\approx 0.7\mathrm{m}$) which is not very well seen in the PIV results.

Figure 11 shows that similar flow pattern was observed also with URANS calculation with the SST model. In the next text, we focus the attention only on the comparison of LBM with the results of SAS model and with the experimental data, and we refer readers to [2] for detailed comparison of the SAS and SST models.

Figure 12, Figure 13 and Figure 14 show the flow field at the vertical cut at $75\mathrm{mm}$ from the symmetry plane. One can again recognize the dominant vortex near the bottom at $x\approx 0.35\mathrm{m}$ and another one at the upper side of the jet at $x\approx 0.7\mathrm{m}$.

The comparison with Figure 15, Figure 16 and Figure 17 corresponding to cut at $z=-75\mathrm{m}\mathrm{m}$ reveals a slight nonsymmetry of the flow field. The vortex at $x\approx 0.05\mathrm{m}$ completely disappears in the FVM simulation and LBM shifts the position of stagnation point a bit upstream with respect to cut at $z=0.75\mathrm{m}$.

Figure 18 and Figure 19 show the flow field at $z=141\mathrm{m}\mathrm{m}$, i.e., close to side wall of the test channel. FVM predicts a separation line at $x\approx 0.2\mathrm{m}$ which is not seen in the results of LBM. Both numerical methods predict a singularity near the bottom at $x\approx 0.5\mathrm{m}$.

The separation line is not visible in the cut at $z=-141\mathrm{m}\mathrm{m}$ where both methods show a vortex in the left part of the test channel, see Figure 20 and Figure 21.

The nonsymmetry is much better visible in the horizontal cuts at $10\mathrm{m}\mathrm{m}$ above the bottom of the channel (Figure 22, Figure 23 and Figure 24). The position of the critical point is shifted a bit from the centerline both for the LBM results and PIV data. Please note that the position of this nodal point is unstable and therefore its interpretation from the time-averaged data is rather inappropriate.

Figure 25, Figure 26 and Figure 27 show the solution in the horizontal cut at $135\mathrm{mm}$ above the channel bottom. One can again see a slight nonsymmetry of the solution both in the velocity magnitude and the vector lines. The FVM method predicts a horizontal vortex in the upstream part of the channel at $z\approx 0.05\mathrm{m}$ whereas LBM shows the vortex at $z\approx -0.05\mathrm{m}$.

Figure 28, Figure 29 and Figure 30 show the structure of the flow field at $y=170\mathrm{mm}$. One can see two non-symmetric vortices at $x\approx 0.45\mathrm{m}$ and $z\approx \pm 0.1\mathrm{m}$ in the LBM results. Similar non-symmetric structure is seen also in FVM calculation and PIV data, although the position is slightly different. The nonsymmetry is also visible in the front part of the channel in the LBM and FVM data.

Very important quantity for engineers is the total pressure loss. The pressure (in lattice units) can be obtained from the density in LBM using the equation of state

$$p=c_s^2\rho ,$$

(34)

where $c_s^2=1/3$. The total pressure is then

$$p_{tot}=p+\frac{1}{2}\rho \left|\right|\overrightarrow{u}{\left|\right|}^2+\rho gy.$$

(35)

The total pressure loss is then calculated from two area averaged values at $x_1=-0.6\mathrm{m}$ and $x_2=2.5\mathrm{m}$ as

$$\Delta p_{tot}=p_{tot,1}-p_{tot,2}.$$

(36)

The LBM method gives the total pressure loss $\Delta p_{tot}=1086.82\mathrm{Pa}$. The same can be evaluated also from FVM which gives $\Delta p_{tot}=1121.71\mathrm{Pa}$. The relative error between LBM and FVM is less than 3.1%. Unfortunately, the experimental device was not instrumented with pressure probes and we cannot therefore compare the numerical data with experimental one.

4. Discussion

The goal of this paper was to investigate the usability of free-surface LBM for the simulation of flows in a model discharge object with the siphon. The attention was focused namely to the structure of time-averaged flow field. The results obtained with the free-surface variant of LBM show very good agreement both with the experimental data obtained with PIV method and with the numerical simulations obtained with the conventional finite volume method published in [2]. Please note that the LBM simulations used several crude simplifications, namely the staircase approximation of the geometry without any near wall mesh refinement or any wall functions and the single relaxation time approach combined with the Smagorinsky LES model.

The prominent feature of LBM is the Cartesian mesh approach which greatly simplifies the preprocessing phase of simulations by removing the lengthy mesh generation process. Moreover, LBM is suitable for implementation on modern parallel computers and GPUs; however, the current LBM simulations have been run on a multicore parallel computer with 72 CPU cores. The required wall time for 1 $\mathrm{s}$ of physical time was approximately $5.5$ $\mathrm{h}$ with almost 80 × 10⁶ cells. The FVM required about 17 $\mathrm{h}$ using 32 cores with approximately 10 × 10⁶ cells. However, it is difficult to directly compare the performance of LBM and FVM because of different mesh size (LBM has almost 8x more elements) or local element sizes (FVM uses mesh refinement near the wall whereas LBM uses a uniform mesh).

The current approach does not contain any bubble model and therefore it is not usable for the simulation of siphon fill-in which has been carried out with FVM in [3]. It can be clearly seen in Figure 4 and Figure 5 where some non-physical “bubbles” remain in the siphon. Next, the lack of a carefully validated wall model in the current approach probably impairs the prediction of the wall shear stress and consequently affects the structure of the flow field near the walls.

However the free-surface LBM was already validated by many authors, the validation was done mostly for relatively simple geometry (dam break, falling droplet, wave breaking). The current study shows the validation for flows in siphon. Moreover, it also shows the importance of proper initialization or a necessity of more advanced model for simulations of transient siphon fill-in.

A more advanced model including either the bubble model such as [28] or a two-fluid model should be investigated in the future. Similarly, the use of regularized LBM or multiple relaxation time variants of LBM should be investigated in combination with more advanced LES models such as WALE model [32].