A Smoothed Particle Hydrodynamics Approach for One-Dimensional Dam Break Flow Simulation with Boussinesq Equations

The critical parameters in a DBF analysis are the surge wave's depth and the wave's arrival time at a downstream location. Traditionally, DBF has been studied using a numerical solution of the Saint Venant equations, which considers hydrostatic pressure distribution along the water depth. However, the distribution of pressure immediately after the dam's failure is non-hydrostatic [1], [2].

Several open channel flows have been studied by numerical solutions of Boussinesq equations using the Eulerian approach [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16]. DBF was simulated numerically by Mohapatra and Chaudhry [6] using Boussinesq equations. Frazao and Guinot [7] developed a hybrid scheme for Boussinesq-type waves in rectangular channels. Shirai et al. [8] investigated the Boussinesq models for their applicability to simulate water wave propagation. Wang et al. [11] compared shallow water and Boussinesq models for cascading DBFs. Magdalena et al. [12] applied the Boussinesq model to simulate DBF and compared the results with the finite volume solution. Darvishi et al. [13] solved the Boussinesq equations for flows over steps and structures. Devkota and Imberger [16] have used the fractional step method to solve Boussinesq equations. Aureli et al. [17] provided the review of DBF modeling for their applications in hydraulic engineering. Some open channel flow studies by Boussinesq equations have also been reported using the Lagrangian approach. Lagrangian method, especially the SPH technique, being a particle-based discretization scheme, is more suitable in dealing with flows with large deformation, such as breaking waves and DBFs. It can efficiently simulate moving boundaries and free surface flows. Since DBF is characterized by large deformations, the SPH method is selected to simulate various dam break scenarios. Chang et al. [18] have used the SPH method to solve Boussinesq equations. However, they considered the first two of the three Boussinesq terms. None of the studies in the literature found that solves Boussinesq equations considering all three terms utilizing the SPH framework.

Researchers [7], [19], [20], [21], [22], [23], [24], [25], [26], [27] have used the SPH method to solve Saint Venant equations. Frazao and Guinot [7] developed a parallel SPH scheme for free surface flows. Researchers [17], [18], [23], [28] presented an approach to model shallow water flows in open channels using SPH. Wang and Shen [19] simulated 1D shallow water DBF using SPH. Ata and Soulaïmani [21] developed a stabilized SPH method for inviscid shallow water flows. De-Leffe et al. [22] modeled shallow water coastal flows by the SPH method. Kao and Chang [24] modeled DBF and flood inundation using the SPH technique. Lin et al. [25] developed a hybrid SPH-Boussinesq model to predict the lifecycle of landslide-generated waves. Researchers [20], [26], [27] developed advanced algorithms to simulate shallow water open channel flows using SPH methodology. Recently, Diwakar and Mohapatra [28] simulated 1D steady and unsteady open channel flows using the SPH model.

This paper numerically solves one-dimensional Saint Venant and Boussinesq equations using the SPH method to simulate 1D DBFs. The present study aims to determine the predictive capability of the SPH model for the effects of Boussinesq terms on the solutions to DBF problems. Therefore, numerical solutions are obtained and compared for the models with and without Boussinesq terms. The predictive capability of the model is obtained by incorporating the Boussinesq terms into the SPH model individually and all together. This study considers all three terms of the Boussinesq equations, and the effects of individual Boussinesq terms for DBF simulations are quantified. Such a study on the SPH methodology considering all the Boussinesq terms is not reported in the literature.

One-dimensional Boussinesq equations for shallow water flows are [29]:

where, $x$=longitudinal direction, $t$=time, $h$=flow depth, $u$=depth-averaged velocity, $g$=acceleration due to gravity, $S_o$=bed slope and $S_f$=friction slope. $B_1$, $B_2$, and $B_3$ are the Boussinesq terms and are given by [29]:

The inclusion of the Boussinesq terms in the flow equations allows for non-hydrostatic pressure distribution. It is introduced by including the second-order derivative of pressure distribution along the water depth. The Saint Venant equations can be assumed to be a particular case of Boussinesq equations by considering Boussinesq terms equal to zero. For detailed explanations and derivations of the Boussinesq equations, readers may refer to the study [29].

The SPH method with an explicit time integration scheme is used to solve the governing equations. First, an intermediate flow field is determined considering the Boussinesq terms, $B_2$, and $B_3$ only. Then, the intermediate flow velocity is corrected using Eq. (4). The solutions are obtained using SPH methodology using self-written code on MATLAB software.

SPH Methodology and Implementation: Following SPH formulation of the study [19], Boussinesq equations for particle $i$ excluding the term $B_1$ can be written as

and,

where, $N$ is the number of particles within the support domain $l$, contributing to the summation having positions $r_j$, $V$ is the volume of the particles, $W$ is the kernel function, and $D/Dt$ refers to the total derivative.

The particle masses are conserved as SPH follows the Lagrangian kinematic approach. Thus, the continuity equation is implicitly satisfied. Therefore, water depth $h$ can also be computed using an SPH approximation (Eq. (7)).

Boussinesq terms, $B_2$, and $B_3$, for particle $i$ are approximated as

$B_1$ is approximated as

where, $\tilde{\mathrm{u}}$ represents the intermediate velocity, and $\mathrm{u}$ denotes the values at the known time level.

The corrected velocity is obtained from the intermediate velocity by

Kernel Function: The following cubic spline function is used as the kernel function.

where,

$$q=\frac{\left\|\mathrm{r}_i-\mathrm{r}_j\right\|}{l}$$

The cubic spline function is accurate and efficient compared to other kernel functions [30].

Artificial Viscosity: Artificial viscosity, $\Pi$, as proposed by the study [31], is used to avoid the interpenetration of particles and to dampen the oscillations.

where, $\alpha$ and $\beta$ are constants, $\bar{c}_{i j}$ is the average sound speed associated with particles $i$ and $j$ and $\mu_{i j}=$ $l\left(\mathrm{u}_i-u_j\right) \cdot\left(\mathrm{r}_i-\mathrm{r}_j\right) /\left[\left(\mathrm{r}_i-\mathrm{r}_j\right)^2+\varepsilon^2\right]$. The term involving $\alpha$ introduces shear and bulk viscosity, and the term involving $\beta$ handles the shock. $\varepsilon$ is a small numerical constant to avoid division by zero.

Source Term: The bed slope and the friction slope terms in Eq. (2) are modeled by Eqs. (14) and (15), respectively.

where, $n$ is the Manning's roughness coefficient, and $b$ represents the bed elevation.

Support Length: Each particle has its support length, which can vary in space and time [19].

in which, $l_0$ and $h_0$ are the initial support length and water depth, respectively.

Boundary Conditions: The SPH particles do not remain in a fixed position due to the Lagrangian nature. This study uses virtual or ghost-type boundary particles [32] to enforce the boundary conditions. Boundary particles are added outside the boundary up to a distance equal to the support radius of the fluid particles. The position of these boundary particles is kept fixed and used only in the summation for particle approximation. The values of the tangential component of the velocity of a virtual particle are taken equal to that of the nearest fluid particle, and the normal component of the velocity is taken opposite to that of the nearest fluid particle. Other properties associated with these boundary particles are taken as equal to the magnitude of the same property of the nearest fluid particle.

Time Integration: An explicit leap-frog time discretization technique integrates particle position and velocity with time. In this scheme, the velocity and the position of a particle $i$ can be obtained as follows:

where, $\Delta t$ is the time step, $\mathrm{r}$ is the position, $\mathrm{u}$ is the velocity at the known time level, and $a$ is the acceleration.

The time step, $\Delta t$, is computed by satisfying the Courant–Friedrichs–Lewy (CFL) stability condition:

where, $C_N$ is the courant number and $c$ refers to the sound speed.

The presented numerical model is applied to simulate various hypothetical DBFs. Input parameters used for the purpose are the length of the channel, $L$=200m; the location of the dam is at 100m from the upstream end of the reservoir; and the initial water depth in the reservoir, $h_u$=1.0m. Different depth ratios [33], [34], $r$, are used in the simulations. The channel bed is considered smooth, and the bed is taken horizontally in these simulations. The flow is simulated using the initial particle number (i.e., number of particles used in the simulation) equal to 4000, and the time step is determined by using Eq. (26) and $C_N$=0.5. All simulations are performed with and without the Boussinesq terms. The initial particle number is chosen based on the sensitivity and convergence analysis (reported elsewhere) and $C_N$=0.5 is considered to have a stable solution.

Four distinct initial water depth ratios (i.e., $r$=0.1, 0.3, 0.5, and 0.7) are considered. The simulated water surface profiles at $t$=15s show that all four zones are similar to those in Saint Venant Equations (Figure 1). However, a lower value of the constant flow depth (the zone before the positive shock front) is obtained in the presence of the Boussinesq terms. In addition, the distance traveled by the positive shock front is higher. For example, the decrease in the flow depth and the increase in distance traveled are 7 percent and 5 percent, respectively, for $r$=0.5. The difference increases as the depth ratios increase. In addition, there are oscillations in the constant flow depth zone when the depth ratio, $r \geq 0.4$.

Evolution of the surface profiles for different depth ratios is assessed for $t$=0 to 2.0s (Figure 2, Figure 3, Figure 4, Figure 5). There is only a marginal effect of the Boussinesq terms on the surface profile in the initial phases ($t \leq 1$s) when $r < 0.4$ (Figure 2 and Figure 3). However, the effect is indicated in the surface profile as time progresses. As shown in Figure 4 and Figure 5, the effects of the Boussinesq terms for higher depth ratios ($r \geq 0.4$) are prominent in the simulations even from the beginning of the DBF.

The individual contribution of the Boussinesq terms on the surface profile is quantified based on the difference in flow depths obtained by including and excluding the Boussinesq terms. Results at two-time instants for $r$=0.5 are presented in Figure 6. Two different time instants are considered for the purpose. The first Boussinesq term is essential (15 percent at $t$=0.5s) at the beginning of the flow. However, the third Boussinesq term is the most influential as time progresses (20 percent at $t$=2s).

A sensitivity analysis concerning the effects of artificial viscosity, bed roughness, and bed slope on surface profiles indicates that (i) as $\alpha$ is increased, the fluctuations in the surface profile are smoothened; (ii) $\beta$ has nominal effects; (iii) there is a retardation of surge waves with an increase in $n$, and (iv) there is a marginal effect of bed slope on the downstream surge propagation.

An explicit numerical procedure using the SPH method is employed to solve the Boussinesq equations for simulating DBF. The following are the main conclusions from the present study:

(1) The DBF surface profile's constant water depth region has a smaller water depth than that obtained by Saint Venant equations.

(2) The positive flood wave propagates faster in the presence of Boussinesq terms.

(3) The effects of Boussinesq terms are insignificant for smaller depth ratios ($r<$ 0.4) in the initial phases. However, it is always prominent for higher depth ratios ($r \geq 0.4$).

(4) Out of the three Boussinesq terms, $B_1$ and $B_3$ have more effects initially. However, $B_3$ influences the results as time progresses. $B_2$ has only a marginal effect.

The presented research can be extended in the following areas:

(1) The present study shows the applicability of the Boussinesq model to DBF scenarios using SPH methodology. The research can further be extended to other flow problems, such as wave breaking and hydraulic jumps where non-hydrostatic pressure distribution is predominant.

(2) The Boussinesq-SPH model is analyzed for 1D flow conditions in the present work. The study on two-dimensional (2D) flows can be performed using a similar approach.

(3) A more detailed and comprehensive study of the presented approach can be performed by evaluating the model performance based on various statistical parameters.