Abstract
The numerical simulation of fast-moving fronts originating from dam or levee breaches is a challenging task for small scale engineering projects. In this work, the use of fully three-dimensional Navier-Stokes (NS) equations and lattice Boltzmann method (LBM) is proposed for testing the validity of, respectively, macroscopic and mesoscopic mathematical models. Macroscopic simulations are performed employing an open-source computational fluid dynamics (CFD) code that solves the NS combined with the volume of fluid (VOF) multiphase method to represent free-surface flows. The mesoscopic model is a front-tracking experimental variant of the LBM. In the proposed LBM the air-gas interface is represented as a surface with zero thickness that handles the passage of the density field from the light to the dense phase and vice versa. A single set of LBM equations represents the liquid phase, while the free surface is characterized by an additional variable, the liquid volume fraction. Case studies show advantages and disadvantages of the proposed LBM and NS with specific regard to the computational efficiency and accuracy in dealing with the simulation of flows through complex geometries. In particular, the validation of the model application is developed by simulating the flow propagating through a synthetic urban setting and comparing results with analytical and experimental laboratory measurements.
1. Introduction
The detailed three-dimensional (3D) simulation of free-surface flows is a challenging task for engineering projects [1] dealing with open flow hydraulics at both the large and small scale because of the computational burden and the large parameterization needed for representing the physics, geometry, and boundary conditions. Nevertheless, in recent times, the always larger availability of remotely sensed data, providing digital topographic and land cover information, and the decreasing impact of the computational burden make the application of small-scale high-resolution 3D numerical models feasible and attractive for engineering studies. In particular, it is foreseen that 3D hydraulic simulations shall support not only fluid mechanics analyses for engineering design at the small scale, but also large-scale applications in urban planning, river restoration, and flood mitigation and management projects.
There are several recent works in the literature that report the ability of the Navier-Stokes (NS) method using the two-dimensional (2D) shallow water (SW) equations [2] to simulate river surface flows, even if the SW hypothesis is not always valid for the inertial effects and significant gradients [3–5]. In fact, several fully three-dimensional NS simulations are available in the literature for representing not only dam-breaks [6, 7], but also unconfined surface flows along straight river segments, bends and confluences [8], as well as impulse waves and tsunamis [4, 9], and water falls [10] among others.
The lattice Boltzmann method (LBM) is an alternative numerical fluid dynamics scheme based on Boltzmann’s kinetic equation [11] that represents the flow dynamics at the macroscopic level by incorporating a microscopic kinetic approach that also preserves the conservation law [12, 13]. LBM has been adopted for the characterization of flows in porous media [14] and for multiphase flows [15–17] among others [18, 19]. The LBM demonstrated to work efficiently for describing the complex physics of non ideal flows and complex geometries [20]. There are also several works investigating the performances of the LBM in simulating shallow free-surface flows [21, 22], but with no emphasis on their use for practical engineering applications.
In this work a fully 3D front-tracking formulation of the LBM, that was originally presented [23] for foam dynamics and that was revisited for free-surface flow applications [24–26], is implemented, validated, and compared to the standard NS approach for transient hydraulic flows, including a synthetic case of flow propagation through an urban setting.
The storage effect of urban areas as well as the hydrograph separation effect resulting in different water transmission rates was investigated in several works [27]. Experiments of dam breaks were conducted [28] to investigate the effect of a flood flow propagating to inundate a single building [29]. Nevertheless, the application and potential validation of LBM for urban modeling of flash floods [30] has not been adequately developed yet.
The novel aspect of the proposed work is the investigation for the potential of the LBM to simulate the flooding in highly urbanized areas with specific regard to the detailed characterization of instantaneous discharge variations, like the one caused by dam or dyke breaks or by intense rainfall of flash floods, producing multidimensional street surface flow paths, which is one of the most complex and difficult hazardous phenomena to manage, prevent, and control. The interaction of water flow with anthropogenic features produces complex 3D flow structures that shall be simulated by taking into consideration all the resulting flow physics such as hydraulics jumps, flow constrictions, and back pressure effects. Disregarding those geometric and hydraulic features would result in a misinterpretation of the estimated water levels providing an inaccurate flood hazard and risk characterization.
In this study, given the impossibility of gathering measured data on such phenomena, a comparative analysis of the performance of LBM and NS in simulating flash flood propagation in urban areas is investigated by a comparison with experimental data of dam break hydraulic effects on a synthetic case study represented by an ideal city.
The present document is organized as follows: in the next section, the theoretical and numerical specifications of both LBM and NS are introduced and described. Then, in the case study section, the results of the application of the proposed methods are inserted with specific regard to the test case of the LBM results compared to an analytical solution; the analysis and comparison of LBM and NS codes applied on the simulation of a rapidly varying flow along a simple geometry of a straight river reach; the evaluation and validation of the LBM and NS code application on the flow propagating through a synthetic urban domain by means of comparison with experimental data of a physical model. The last two sections represent, respectively, the discussion and conclusive remarks.
2. Methods
2.1. The Navier-Stokes Numerical Scheme
The most general macroscopic model for hydraulic flows may be represented by the incompressible Navier-Stokes equations as follows: where is the fluid density, is a specific prescribed body force (i.e., gravity force), and is the fluid dynamic viscosity.
As in hydraulic flows, at least two phases are always present, water and air; the above system of equations must be applied to a multiphase system. The simulation of multiphase flows is a challenging task due to the inherent complexity of the involved phenomena (i.e., moving interfaces with complex topology) and represents one of the leading edges of computational physics.
The following two main approaches are widely used to simulate multiphase flows: the Euler-Lagrange approach and the Euler-Euler approach. The latter approach, usually chosen in hydraulics applications, is based on the assumption that the volume of a phase cannot be occupied by other phases, thus introducing the concept of phase volume fractions as continuous functions of space and time.
For the present work, we employed the so-called volume of fluid (VOF) method, [31, 32] designed for two or more immiscible fluids where only one fluid (i.e., air) is compressible and the position of the interface between the fluids is of interest. The VOF method is a surface-tracking technique applied to a fixed Eulerian mesh, in which a specie transport equation is used to determine the relative volume fraction of the two phases, or phase fraction, in each computational cell. Practically, a single set of Reynolds-averaged Navier-Stokes equations is solved and shared by the fluids, and, for the additional phase, its volume fraction is tracked throughout the domain.
Therefore, the full set of governing equations for the fluid flow are as follows: where is the velocity vector field, is the pressure field, is the turbulent eddy viscosity, is the strain rate tensor defined by , is the surface tension, and is the surface curvature.
The nature of the VOF method means that an interface between the species is not explicitly computed, but rather emerges as a property of the phase fraction field. Since the phase fraction can have any value between 0 and 1, the interface is never sharply defined, but occupies a volume around the region where a sharp interface should exist.
The model used in this work is based on an open source computational fluid dynamics (CFD) platform named OpenFOAM [33], freely available on the Internet. OpenFOAM, primarily designed for problems in continuum mechanics, uses the tensorial approach and object-oriented techniques. It provides a fundamental platform to write new solvers for different problems as long as the problem can be written in tensorial partial differential equation form.
2.2. The LBM Numerical Scheme
This section introduces the theoretical background of the proposed LBM with specific regard to the following:(i) the single-phase LBM and the boundary conditions,(ii) the front tracking concept,that are, respectively, inserted in the two following sections.
2.2.1. Single-Phase Lattice Boltzmann Method and Boundary Conditions
Artificial and natural fluvial hydrodynamics may be simulated at the macroscopic scale using the LBM that is structured as a mesoscopic (i.e., scaled larger than micro and smaller than macro) method. LBM treats the water mass by analysing its density probability function, , where is the particle spatial location within the lattice domain, is the time variable, and is the kinematic component. The dynamic of the particle distribution function is described in (3) that represents a discrete lattice Boltzmann equation, representing the evolution in time of the fluid flow by providing a numerical solution of the modification of the density probability function: where the left and right terms represent, respectively, the molecular free-streaming and the particle collisions. is an additional term accounting for internal (e.g., collision) and external forcing processes (e.g., electric and magnetic fields or gravity).
For water-driven processes, it is reasonable to assume a discrete set of directions along which particles may route and, as a consequence, is restricted to a maximum of 7 to 9 for 2D or 19 directions for 3D flows. This 19-direction discretized cubic lattice model implemented here (namely, D3Q19) is characterized by the set of predefined velocity vectors schematically represented in Figure 1 and reported in Table 1.
The collision component is integrated by implementing the Bhatnagar-Gross-Krook Model (BGK) [34] that characterizes the density modifications caused by particle collisions as related to dynamic modifications from the local equilibrium [35, 36].
If the relative Knudsen and Mach numbers are small enough, this formulation is able to represent the dynamics of a fluid with pressure and a kinematic viscosity , where is the lattice sound velocity [12].
The proposed numerical approach is implemented on a structured Cartesian lattice (i.e., mesh) where the particles may only reside on the mesh nodes.
In the proposed model boundary nodes are solved using the bounce-back scheme meaning that particles colliding with a stationary wall simply reverse their momentum [12, 37]. In particular, the half-way bounce-back scheme is here applied by posing the reflecting wall at a medium distance between boundary nodes and the physical walls [38]. As a result, the proposed model can simulate both solid features (e.g., walls, levees, etc.) of known geometry or open boundaries of predefined velocity, discharge, or pressure function.
The turbulence flow field is implemented in the proposed LBM model by employing a Large Eddy Simulation (LES) approach [39] that characterizes the flow field simulating large-scale structures (resolved grid scales), which are assumed to be mostly predominating in defining the transport of mass and momentum, and not directly simulating the small-scale local effects. The definition of the boundary between small and large processes is based on the Smagorinsky model [40].
2.2.2. Front-Tracking Model
In the Front-Tracking (FT) formulation proposed here [23], the interface is treated as a zero thickness, mathematical surface, across which the density field jumps from the light to the dense phase and vice versa and is tracked throughout the computational domain through an additional scalar variable (i.e., liquid mass fraction). A single set of LB equations is solved for both the heavy fluids with respect to the light fluids (i.e., gas) in agreement with the front-tracking methodology that do not aim at providing any description of the physics of the phase-transition between the two phases. The main difference, which also represents the main advantage of the LBM front-tracking approach, with respect to the traditional continuum methods, is that the single-species transport equation, used in Navier-Stokes/VOF model to determine the relative volume fraction of the two phases in each computational cell, is not needed, since the free-surface tracking is automatically performed by advancing the fluid elements. In practical terms, only the liquid phase flow is numerically solved, while the gas characterization is neglected, and an additional variable (i.e., the liquid volume fraction) is structured to track the free surface throughout the computational domain. This variable is numerically solved using the following equation: in which , and represent respectively empty (gas), liquid, and interface cells.
The computation scheme is integrated numerically in three steps as follows: interface positioning and flow mass motion update defined by analysing cell fluxes; definition of the boundary conditions at the interface to separate the gas and liquid state; update of cell type. As a result, interface cells defined by the position and time are defined as follows: in which characterizes the postcollision state and the sign is related to the mass direction where has opposite direction with respect to . Mass dynamics in space and time are defined as
Combining (6) and (7), the continuity principium is in respect with the condition that direct modification of state, from liquid to gas or vice versa, is not permitted outside the interface. In fact, liquid and gas cells are only allowed to transform into interface cells, whereas interface cells can be transformed into both gas and liquid cells. In this way, mass and density are completely decoupled and, as a results, mass evolution does not affect the particle distribution functions.
3. Numerical Results for Three Case Studies
The proposed LBM and NS are applied for the following three test cases:(a) the free-surface flow routed along a straight channel, with LBM modeling results that are validated by comparison with two different types of analytical solutions;(b) the free-surface flow routed along a straight channel deriving from the partial asymmetrical collapse of a vertical dam embankment with results that are compared in terms of hydraulics and computational efficiency;(c) the simulation of a severe transient flow through a synthetic urban district.The three tests are hereafter presented in detail.
3.1. Test (a): Fast Flow Propagating along a Straight Channel
The routing of a fast moving flow resulting from a dam break, propagated through a 20 m long straight channel (see Figure 2), is simulated, and the numerical results are compared to the analytical solution provided by Ritter [41].
The channel bed is horizontal. Bottom and wall friction is set to null. The numerical mesh is made up of 2,000 cubic elements with a grid-spacing of 0.01 m. At time zero, a 0.5 m by 7 m volume of water ideally positioned at the most, upstream section of the channel is released instantaneously and two resulting waves propagate, respectively, downstream and upstream back to the reservoir. At a given time from the start of the simulation (i.e., dam break/breach or ideal vertical wall removal), three different regions are defined (Figure 2) as follows:(1)for , we have m,(2)for , the water surface profile is parabolic (3)for we have ,where is the water elevation. Figure 3 shows the agreement between numerical and analytical data apart some minor expected differences at the boundaries [42]. In fact, where Ritter's analytical scheme is not suited for this comparison, other methods were selected like the Chanson’s solution [43]: that is also plotted in Figure 3.
To further prove the validity of the code, a dam break flow propagating along a 40 m long straight channel is simulated. At time zero, a 0.23 m by 18 m water volume is released instantaneously. In Figure 4, LBM numerical results are compared to the Ritter [41] and Chanson [43] as well as to the experimental measurements provided by [44]. The comparison is performed in terms of dimensionless instantaneous free-surface profiles between numerical and analytical data (assuming ).
(a)
(b)
3.2. Test (b): Asymmetrical Dam Break—LBM and NS Numerical Models
The second test is characterized by the simulation of a submersion wave caused by the partial collapse of a nonsymmetrical sluice gate (or partial asymmetric dam breach). The spatial 3D computational domain is where the boundaries of the domain are impermeable walls. The nonuniform initial condition for the liquid phase fraction is specified: the water surface level is initially 10 m for the upstream section and 5 m for the downstream section. An unsteady flow is generated by the instantaneous collapse of an asymmetrical 75 m long portion of the dam (see Figure 5). The bottom is flat, and ground resistance to the motion is neglected.
Numerical results of several authors are available in the literature for this case [6, 45, 46], that is, a standard benchmark for hydraulic models in dam break applications. We chose to compare the LBM and NS, described in the previous section, in terms of computational efficiency. Calculations have been performed with a time step of 0.01 s on a computational mesh of 800,000 nodes (i.e., lattice sites for the LBM). Figure 6 shows the comparison of the computed water level 7.2 seconds after the breach, when the flow reaches the left side of the tank, while Figure 7 shows the simulated water surface levels along sections A-A (longitudinal m) and B-B (transverse m). Water level hydrographs at selected points of the domain have been also evaluated, and Figure 8 shows the results obtained for points P1 and P3 with both LBM and NS approaches.
(a)
(b)
(a)
(b)
(a)
(b)
Results show that NS and LBM, using the same grid resolution, are in good agreement in terms of water heights and hydrograph, along the section AA, while the results for the section BB (Figure 7(b)) are remarkably different especially in terms of water depths. This is expected with specific regard to the central part of the domain where the dam break originates and the flow is significantly 3D, while in the other portions of the domain, through which the emptying of the system applies (from the noncollapsing area towards where the breach is) in the direction of the flow that is parallel to the gate, the two models give different results with greater water depths evaluated by NS.
From the computational point of view, using the same hardware, the LBM model has completed the tests in one-fifth of the time needed by the NS. Results are summarized in Table 2 that also shows that the LBM code is able to provide a higher level of details of the flow pattern using the same computational time or to provide the same accuracy and details in less time. Figure 9, synthesizing the LBM versus NS computational efficiency comparison test, shows that the LBM model can use a mesh 65% denser than NS one ( cm, m) for the same computational time.
3.3. Test (c): Severe Transient Flow through a Synthetic Urban District
This test presents the simulation of the effects of a severe transient flow on a synthetic urban area. This case study replicates the experiment developed at the Hydraulic Laboratory of the Civil and Environmental Engineering Department of the Université Catholique de Louvain in Belgium [47]. The geometric configuration is characterized by a regular distribution of rectangular obstacles (buildings) intercepting the flow in geometric domain. The realistic average ratio between building and street widths of 1 to 3 was selected by interpreting a real urban configuration (the aerial view of Brussels in Belgium). Buildingsʼ heights are simulated as high enough in order not to be overtopped by the flow. The experimental physical setup consists of a 36 m long by 3.6 m wide channel, with a flat profile (Figure 10). Buildings, wooden parallelepipeds of 30 cm by 30 cm, and streets that are 10 cm wide are positioned along the flow direction. The initial reservoir water level is 0.40 m upstream. Water surface measurements were performed [47] by means of several resistive level gauges for estimating the water surface geometric and dynamic properties (depths and surface velocity fields).
The aim of this test is to investigate the potential of LBM and NS to correctly represent the complex dynamics of transient flows. In particular, the physics of interest is related to those flood phenomena propagated through urban setting for which an initial severe transient phase is often followed by a slow transient that could approximate under certain steady state conditions. Nevertheless, since the first instants correspond to a severe transient flow involving specific flow features characterized by high velocities and transcritical flows, wave impacts on structures, or higher water rise following wave reflection against a structure, the LBM and NS are challenged to provide such unsteady irregular behavior followed by the regular steady conditions of the following sequences.
3.3.1. Simulation Setup—Mesoscopic Approach Using LBM
Due to the large-scale of the domain in horizontal plane, the test case has been replicated using a computational mesh of cubes of 0.025 m resolution producing a total number of 4,976,640 nodes. The downstream region has been represented assigning a liquid volume fraction of 0.5 to the first vertical row of nodes. Considering the dimension of the domain, the mesh is quite rough and only few computational nodes describe the streets that are only 0.1 m wide.
3.3.2. Simulation Setup—Continuous Approach Using NS
For all simulations, the computational domain is the same with the following dimensions: 36 m long, 3.6 m wide, and 1 m high. Five scenarios have been developed for the NS simulation, varying the resolution of the mesh in the three directions. The computational mesh for each scenario derives from the intersection and refinement of the mesh composed of hexahedra (configuration A, see Table 3) with the 3D surface of the synthetic urban setting (Figure 11). In configuration A a grid variation along the vertical with a ratio of 1 to 10 is used. The result of such a refinement process is represented by configuration B (see Table 3). For the five scenarios, M1 to M5, the diverse grid resolutions lead to discrepancies between the numerical results and physical measurements both inside and outside the city. The comparison between experimental data and numerical results provided by the LBM and the NS, scenarios M1, M2, M5, for the section m at = 5, 6 and 10 s after the break are presented in Figure 12.
(a)
(b)
(c)
(a)
(b)
(c)
The free-surface profile includes a hydraulic jump at 5 s after the break. In the graphs for = 5 s and = 6 s, when the flow depth in the city is still low, the flow reflection against the buildings is clearly visible for both LBM and NS simulations, and the LBM model correctly predicts the free surface at low values as the NS model with the finest mesh. At = 10 s, the upstream hydraulic jump is still present, but the water level has significantly increased in the streets. The flow in the streets evolves from supercritical (with a control section near the street entrance) to subcritical (with a control section at the street exit).
The overall agreement between the numerical and experimental results is promising for the use of the mesoscopic numerical model for replicating the complex physics of rapidly varying flows in such a dense urban area, even if, due to the mesh resolution, results of lattice-based simulations are less accurate than the NS because of the mesh refinement (especially for M5). Analyzing the effect of mesh refinement, the obstacle mesh refinement of M3–M5 is able to better reproduce the unsteady hydraulic features and flow processes like hydraulics jumps, while the coarser meshes of M1-M2 are characterized by some errors where the flow behavior is irregular.
For LBM at the time steps = 5 and 6 s, significant differences appear between the two numerical models. For example, the depression and the hydraulic jump at the entrance of the urban district ( = 5.20 m) is absent in the LBM coarse mesh result. One limitation of the used LBM code is that the computational mesh is not refined locally for example in a single area of interest. This disadvantage with respect to the used NS model clearly indicates the utility of employing mesh refinement techniques, as those developed for standard LBM [48].
Figure 13 shows the results at seconds for 5 different meshes, while in Figure 14 snapshots of free-surface evolution at different time steps (0, 2, 4, 6, 8, and 10 s) are inserted. Figure 15 reports the value of velocity magnitude for the same time steps.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(f)
Observations indicated that the flow rises at the city front before entering the streets, after the wave impact, a phenomenon that is similar to the impact against a single obstacle rather than a series of obstacles. A hydraulic jump is represented at the impact section (Figure 14), with the water level that is locally higher with respect to the same case without buildings, while a wake zone is developed immediately downstream of the city.
4. Conclusions
In this study an experimental fully 3D flow modeling framework is implemented to numerically simulate the unsteady irregular flows originating from dam breaks with specific regard to the complex interaction between the geometrical urban features and the associated physical hydraulic processes. The proposed novel technique is based on a front-tracking variant of the lattice Boltzmann method that is here investigated as a valid alternative to the standard numerical Navier-Stokes model in the treatment of instantaneously varying flows interacting with artificial structures. In particular, test cases evaluate the accuracy and the computational efficiency of numerical simulations with respect to analytical solutions and experimental data. The critical analysis of model performances provides the following conclusive remarks.
The proposed LBM is able to predict with a high level of accuracy typical hydraulic phenomena characterizing engineering applications related to the flooding of simple straight channels (test (a)) and asymmetrical front waves in rectangular basins (test (b)). More specifically, as compared to the fully 3D NS model, it is reported the validity and robustness of LBM results that, while implementing the same resolution of the geometric domain, are demonstrated to be more computationally efficient (test (b)).
The last test investigates the diverse hydraulic unsteady processes, like the hydraulic jump, return, and diverged flow, that originate from the passage of the flood wave through a synthetic urban setting (test (c)). The application of the presented 3D numerical schemes shows the major advantages and disadvantages of both LBM and NS schemes confirming the ability of LBM in representing such abruptly changing and irregular flow dynamics. Nevertheless, the mesoscopic LBM approach is significantly bound by the regular structured mesh, while the continuous NS model performs well, also in large domains, especially when using a more accurate grid (M5) that is adaptively refined to more accurately represent obstacles.