## Abstract

Three-dimensional (3D) hydraulic modelling of rapidly varying surface flows is a challenging task for practical engineering applications. One example is represented by the fast-moving fronts originating from dam breaches that proceed downstream through artificial channels. In this work, a fully 3D lattice Boltzmann method (LBM) is tested. The numerical model is a front-tracking variant of the LBM, being the free surface tracked through the liquid volume fraction. Model performances are evaluated simulating the effect of dam-break flows on synthetic settings schematically represented by an artificial domain and comparing results with analytical data and experimental laboratory measurements. Obtained results are promising for the use of LBM for practical engineering applications.

- dam break
- free-surface flows
- hydraulic modelling
- lattice Boltzmann model

## INTRODUCTION

Hydraulic models for simulating free-surface flows using three-dimensional (3D) schemes have been increasingly used in the last decade for several reasons. The large availability of high-resolution 3D digital topographic information on the one hand and the increasing computer power and efficiency on the other have paved the way for scientists to develop and test always more complex modelling algorithms. Moreover, the need for detailed deterministic hydraulic simulations for urban and land planning projects (Wu & Chau 2006; Manciola *et al*. 2009) with specific regard to the physically based characterization of the impact of territorial transformations and man-made interventions on the hydraulics of the domain – that means by providing a detailed numerical characterization of the water surface geometry and associated flow depths and velocities – is challenging for hydraulic modellers to accurately simulate the actual conditions with respect to the design (post-intervention) scenario.

At present, standard models for characterizing complex systems are based on the Navier–Stokes continuum equations (NSE). The shallow water equations (SWE) have proven to be computationally efficient and suitable for representing a wide range of free-surface flows (Lai & Khan 2012; Petaccia *et al*. 2013; Pu *et al*. 2013). Nevertheless, flows are often fully 3D with strong inertial effects and high free-surface gradients that break the main hypothesis of the SWE (Chau & Jiang 2001; Stoesser *et al*. 2003; Chau & Jiang 2004; Biscarini 2010; Biscarini *et al*. 2010; Ferrari *et al*. 2010; Filonovich *et al*. 2013).

A possible remarkable alternative to the NSE is represented by the lattice Boltzmann method (LBM), which is based on a discretized Boltzmann kinetic equation (Higuera *et al*. 1989; Qian *et al*. 1992; Ubertini *et al*. 2010). The LBM has been recently adopted for different applications: the simulation of flows in porous media (Ginzburg & D'Humières 2007; Zarghami *et al*. 2014), and for multicomponent and multiphase flows (Falcucci *et al*. 2010) among others. The LBM has been also implemented to simulate free-surface flows by employing the shallow water hypothesis (Zhou 2004; Budinski 2012).

In this technical note, a fully 3D front-tracking (FT) formulation of the LBM is proposed. The method, originally presented by Körner *et al*. (2005) for foam simulations, was further developed by some authors for free-surface hydraulic flows (Biscarini *et al*. 2011, 2013; Falcucci *et al*. 2011). The numerical model is here tested against experimental measurements of a partial symmetrical collapse of a vertical dam embankment. The novelty aspect of the present work with respect to previous findings is mainly represented by the detailed comparison of the transient flow field and free-surface evolution, given the availability of literature experimental measurements at different points of the domain and at different times.

## THE FT LBM NUMERICAL SCHEME

The LBM treats the fluid through its discrete nature, by analysing the dynamics of the density probability functions, , which represent the probability of finding a particle moving with speed **c _{i}**, in the spatial location

**x**, at time

*t*. The particle velocities being restricted to a limited number of directions

*b*, the discrete lattice Boltzmann equation reads as follows:1where the right side represents molecular collisions implemented as a single-time, τ, relaxation to local equilibrium (Bhatnagar

*et al*. 1954) and

**F**accounts for external forces (e.g., gravity).

_{i}The set of discrete speeds must be chosen in order to guarantee mass, momentum and energy conservation, as well as rotational symmetry, in order to recover the correct fluid dynamic equations in the macroscopic limit. Among the possible options proposed for 3D flows, a 19-direction (i.e., *b* = 19) cubic lattice model is implemented here (namely D3Q19), as a good balance between computational reliability and efficiency (Succi 2001). The resulting numerical mesh is a structured cubic Cartesian lattice where the particles may only reside on the mesh nodes.

The local equilibrium function is defined by expanding in Taylor series the Maxwell–Boltzmann distribution2where *w _{i}* are weighting factors, and and are, respectively, the fluid density and velocity. At small Knudsen and Mach numbers, this formulation is able to represent the dynamics of a fluid with pressure and a kinematic viscosity , where is the lattice sound velocity (Succi 2001).

In the proposed model, boundary nodes are solved using the half-way bounce-back scheme (Skordos 1993; Inamuro *et al*. 1995), which allows modelling complex boundaries in terms of elementary mechanical rules (bounce-back, reflections). The bounce-back scheme can be applied to several boundary conditions usually found in hydraulic problems, such as no-sip, free-slip, open flows and moving boundaries. Curvilinear surfaces can be handled through a stair step approximation, which is sufficiently accurate for the typical geometries of the large-scale hydraulic flows discussed herein.

Turbulence is modelled through the large-eddy simulation (LES) approach; in particular, we employ the Smagorinsky–Lilly sub-grid scale with a Smagorinsky constant of 0.15 (Galperin & Orzag 1993).

The free-surface flows are computed through a FT variant of the LBM (Körner *et al*. 2005; Biscarini *et al*. 2013). As the density ratio between water and air is of the order of 10^{3}, the dynamics of the gas phase is neglected, it being assumed that air is constantly at atmospheric pressure, and the interface is treated as a moving zero-thickness boundary across which the density field jumps from the light to the dense phase and vice versa. The free surface is tracked throughout the computational domain by means of the liquid mass fraction3in which *G*, *L* and *I* represent, respectively, empty (gas), liquid and interface cells.

The dynamics of the liquid phase, which governs the flow completely, is modelled through a single set of LB equations, by computing the particle distribution functions of liquid and interface cells, while excluding from the computation the gas state. The liquid mass fraction is determined by mass dynamics in space and time45in which *pc* characterizes the post-collision state and *–i* defines the opposite direction with respect to *i*.

As the gas dynamics is neglected, the distribution functions coming from the gas cells are never accessed and need to be reconstructed. Here, we use the reconstruction scheme proposed by Körner *et al*. (2005) by setting that the water pressure equals the air pressure at the interface.

## RESULTS AND DISCUSSION

Most importantly, a grid convergence analysis is performed by comparing numerical results with analytical findings from a two-dimensional dam break. We consider a 20 m long channel initially filled with water to a height of 0.5 m (*h _{0}*) in the first 7 m of the channel (

*X*

_{0}). The simulation is performed for different grid sizes (0.1, 0.04, 0.02, 0.01, 0.005, 0.0025 m) and the present results are compared with the analytical solution of the problem proposed by Ritter (1892) and Chanson (2009). The water depths (

*h*) of

_{i}*m*I-cells are compared with the analytical solutions (

*h*) and the mean error is calculated as follows:6

_{as}Figure 1(a) shows the error decay with the grid size. Figure 1(b) shows the wave height differences at *t* = 0.5 s simulated by LBM equation in a test point at position *x* = 8 m. The resulting convergence order of the scheme is about 2.9.

The proposed LBM is then applied for the simulation of the transient free-surface flow deriving from the partial symmetrical collapse of a vertical dam embankment. A detailed analysis of the transient flow field is performed by comparing numerical results with laboratory experimental measurements on a dry bed (Fraccarollo & Toro 1995; Ferrari *et al*. 2010). The dam rupture is simulated by a sudden removal of a gate, considered instantaneous. Water in the tank is initially in steady state with a depth equal to 0.6 m. Measurements of pressure, water depth and velocity components for a period of 10 s are available at gauges located in the tank and downstream along the plane bottom at the positions shown in Figure 2. After that time, the water flow slowly varies and reaches the quasi-steady state with a hydrostatic pressure distribution that characterizes almost the entire domain.

Figure 3 shows both numerical and experimental hydrographs at points 4 and 0 along the gate. At the initial stage of the tank emptying (i.e., first 2 seconds), strong three-dimensional inertial effects are clearly visible, as the high free-surface gradient close to the gate produces a local non-negligible acceleration component on the gravity direction. This is more evident at point 0, as point 4 is closer to the wall (i.e., the gate portion that did not collapse).

In Figure 4, the experimental x-components of the velocity, *V _{x}*, are compared with our results for point 0. Eight measurements are provided at different heights from the bottom, starting from 50 mm from the bottom up to 0.4 m, every 50 mm. Both experimental and numerical velocity profiles at point 0, shown in Figure 4, are characterized by a typical shape: a fast (i.e., <1.5 s) quasi-linear increase to the maximum value followed by a monotonic velocity decrease in the longitudinal direction. The ending time of the curves reduces with higher

*z*values, as the water surface elevation decreases after the gate break, leaving the higher gauges outside the flow domain.

## CONCLUSIONS

In this study, a fully 3D FT variant of the LBM is demonstrated to be a valid model in the treatment of instantaneously varying flows interacting with artificial structures. In particular, we demonstrate that the proposed LBM is able to characterize in detail the complex hydraulics of a test case like the symmetrical dam break flows propagating within a rectangular domain. The numerical model is able to predict the experimental hydrographs and velocity profiles at different points of the domain and at different times.

Results are, thus, promising for the use of the proposed LBM for practical engineering applications. Nevertheless, future work should include: the treatment of moving bodies and curvilinear boundaries; the introduction of grid-refinement techniques to further increase computational efficiency; the introduction of additional phases (i.e., solid, landslide material); the development of sensitivity analyses on modelling and geometric parameters (e.g., roughness effects, mesh resolution); the simulation of real case scenarios of unsteady flows through natural topography; and the validation of simulation results with real case measurements.

- First received 18 February 2014.
- Accepted in revised form 8 June 2014.

- © IWA Publishing 2015

Sign-up for alerts