## Abstract

The silted-up sediment in the reservoir may have a significant influence on the propagation of dam-break flows. In this paper, a three-dimensional numerical simulation of the silted-up dam-break flow is carried out. In this paper, simulations of three-dimensional silted-up dam-break flow are carried out. A kind of Eulerian–Eulerian two-fluid model (TFM), coupled level set and volume of fluid (CLSVOF) methods, is presented. In order to calculate the motions of the air–water interface and the sediment simultaneously, kinetic particle theory (KPT) and computational fluid dynamics (CFD) are combined. The rheology-based constitutive equations of sediment are also considered to simulate scouring and deposition. In addition, a partial-slip boundary condition (BC) for the velocity of the sediment phase at stationary walls is implemented. The simulation results of the benchmark case demonstrate that the proposed model can effectively simulate the silted-up dam-break flow while taking into account multi-interface capturing problems. Subsequently, the simulations of the silted-up dam-break flow over dry are investigated numerically in a three-dimensional long channel. The simulated results reveal that, in the dam-break flows, the silted-up sediment height has a significant influence on wave propagation, dynamic pressure loads, sediment transport, and sediment deposition.

## HIGHLIGHTS

Three-dimensional simulation of silted-up dam-break flow is carried out.

The air–water interface movement is captured by CLSVOF method.

A Eulerian–Eulerian multiphase model coupling kinetic particle theory and computational fluid dynamics is used.

The rheology-based constitutive equations of sediment are considered.

The effect of the silted-up sediment height during a dam-break flow is investigated.

## INTRODUCTION

Dam-break flood is triggered by the sudden structural damage of dams, which leads to a surge in downstream water level and generates unstable flood waves. In general, the investigation of dam-break flows could be roughly divided into air–water interface, dynamic pressure loads, and sediment transport. For air–water interface, Kocaman & Ozmen-cagatay (2012) explored the effect of lateral channel contraction on dam-break waves. Kocaman & Ozmen-cagatay (2015) conducted a similar dam-break flow experiment to investigate water level fluctuations and free surface evolution. Ye & Zhao (2017) used the two-liquid volume of fluid (VOF) approach to analyze the evolution of water–water interface between reservoir water and downstream water in wet-bed dam-break flows. By using records of water surface profiles, stage hydrographs, and cross-sectional mean velocities, Wang *et al.* (2020) investigated the similarity of dam-break flows on dry beds. In addition, three-dimensional simulations of dam-break flows striking rigid structures have also been presented by Arote *et al.* (2021) and Alloush *et al.* (2021). These studies enable a thorough understanding of flood front propagation and water stage changes of clear-water dam-break flows.

For dynamic pressure loads in dam-break flows, Lobovský *et al.* (2014) conducted a series of experiments to study the pressure loads on the downstream vertical wall during dam breaks. Furthermore, Mohamed *et al.* (2018) performed an experimental investigation of the dam-break flow impacting on a vertical cylinder positioned over a horizontal dry bed. Impacting on downstream walls or inflexible obstacles would usually lead to violent free surface deformations and surging dynamic pressure loads in dam-break flows (Meng *et al.* 2021a; He *et al.* 2022; Maghsoodi *et al.* 2022). The impacting pressure of dam-break floods on a solitary structure or various types of obstacles in the downstream region has also been numerically studied (Meng *et al.* 2021a; He *et al.* 2022; Maghsoodi *et al.* 2022), which provides accurate fluctuation profiles of air–water interface and time evolution of dynamic pressure on the rigid-body surface.

However, the propagation of air–water–sediment multiphase dam-break flows may be significantly different from that of air–water two-phase flow. The influence of sediment on dam-break wave propagation may be of great importance. In recent years, studies on the scouring process and the meandering effect have also shown that sediment plays an important role in dam-break flows (Bihs & Olsen 2011; Pu & Lim 2014; Zheng *et al.* 2016). In addition, turbulence has been shown to be another key effect to be investigated (Pu *et al.* 2012a, 2012b, 2014; Pu 2015), especially in dam-break flows. Moreover, sediment transport in dam-break flows may result in turbulence and affect the dam-break flow process. Previous studies mainly focused on clear-water dam-break flows. There is still a lack of quantitative investigations of the effect of silted-up sediment height on dam-break-induced flood. However, there has been no numerical investigation on silted-up dam-break flows because it remains a challenge for existing solvers to simulate rheology-based sediment movement with complicated fluctuations in air–water interface. Therefore, to fill the gap and improve numerical simulation methods, this study aims to investigate the effects of silted-up sediment height on a dam-break flow by proposing an Eulerian–Eulerian multiphase model. Different from Lagrangian methods such as SPH (Chang & Chang 2013; Pu *et al.* 2013; Shi *et al.* 2019), this model adopts the Eulerian approach. Compared with the one-fluid Eulerian model such as SWEs (Pu *et al.* 2012a, 2012b), this study employs the two-fluid Eulerian multiphase method to simulate both air–water interface and sediment movement.

The innovations of this study can be summarized as follows:

- (1)
The proposed Eulerian–Eulerian two-fluid model (TFM), coupled with the level set and volume of fluid (CLSVOF) method, is a combination of the kinetic particle theory (KPT) and computational fluid dynamics (CFD) to calculate air–water interface evolution and sediment movement.

- (2)
In the proposed multiphase model, the rheology-based constitutive equations and a partial-slip boundary condition (BC) of sediment are implemented to simulate scouring and deposition.

- (3)
Unlike previous clear-water dam-break flows that only considered the water interface, this study also considers the effect of silted-up sediment, providing a new insight into the emerging silted-up dam-break flow problem.

## EULERIAN–EULERIAN CLSVOF METHOD

### Interface capturing

*et al.*2019). The Multidimensional Universal Limiter with Explicit Solution algorithm in OpenFOAM (Nguyen

*et al.*2020) is employed in the VOF method to counteract ambiguities near the phase interface caused by numerical dissipation:where is the fraction of volume occupied by the water phase and is the fluid velocity vector; is the compression velocity;

*c*is a controllable compression parameter and is set to 1. After solving the above VOF Equation (1), the initial value of the LS function is calculated by , where and is the size of grid cell. Then, the LS function is re-distanced by solving the re-initialization equation (An & Yu 2020):where is the sign function and represents the artificial time step. After solving the above LS Equation (2), the surface tension of the air–water interface can be calculated by , where is the surface tension coefficient, is the interface curvature, and is the Dirac delta function (Yu

*et al.*2019; An & Yu 2020). In addition, the physical properties of the fluid phases (including air and water) will be updated by weighting factors at each time step, calculated by the following formula:where the subscripts

*l*and

*g*represent the water phase and the air phase, respectively; represents density and represents viscosity.

*e*is the diffusion coefficient set at 1.0 × 10

^{−8}. Here, because diffusion enhances numerical stability, Equation (4) is an advection-diffusion equation rather than a pure advection equation, as suggested by Lee (2021).

### Momentum equations

*f*and

*s*represent the fluid phases (including water and air) and the sediment phase, respectively. Here,; ; is the total pressure of the fluid phase; is the gravity vector; is the position vector of the cell center; is the viscous stress of the liquid phase, while is the intergranular stress of the sediment phase; is the sub-particle scale (SPS) stress; and represents the drag force between the sediment phase and the fluid phases (including air and water) (explained in Section 2.4).

*k*

*=*

*f, s*) is the sub-grid scale eddy viscosity (Ferziger & Perić 2002); (

*k*

*=*

*f, s*) represents the rate-of-strain tensor, in which the superscript

*T*represents the transposition of the second-order tensor matrix; and is the unit second-order tensor. The improved Smagorinsky (Ferziger & Perić 2002) model is used to determine the value of (

*k*

*=*

*f, s*):where indicates the length of the filter in Large Eddy Simulation (LES) (Ferziger & Perić 2002) and is defined as the cell size; (

*k*

*=*

*f, s*) is the norm of ;

*N*is a coefficient and is set to be 5; is the maximum sediment volumetric concentration; (

*k*

*=*

*f, s*) represents the model parameter and (

*k*

*=*

*f, s*) is taken.

### Rheology-based constitutive equations of sediment

*et al.*2021b) will influence sediment motion. The intergranular stress is calculated using the equation described by Meng

*et al.*(2021b) and Shi

*et al.*(2019):where exists in the persistent contact, particle inertia, and fluid viscosity behaviors, and represents the elastic effect and is important when sediment is silted in a compact bed.

*et al.*(2021b):where is the jamming volume fraction,

*b*and

*a*are the model parameters whose values are set to be 2.0 and 0.11, respectively, and

*d*is the diameter of the sediment particle.

*K*is the coefficient related to the Young's modulus of the sediment; is the model parameter, and and account for the random loose-packing concentration and the random dense-packing concentration, respectively. As increases to , the component comes into play. When is just slightly less than , this implies the transition from a fluid-like state to a solid-like state.

*I*, and

*m*is a regularization parameter that is set to 50 in this study. The following formula is adopted to compute the friction coefficient (Meng

*et al.*2022):

### Sediment–fluid interaction model

In silted-up dam-break flows, the carrier phases (air and water) are responsible for promoting the motion of the solid phase (sediment), and the sediment will in turn affect the fluid movement as well. Since drag force plays a dominant role in dam-break flows with high sediment concentration (Shi *et al.* 2019; Meng *et al.* 2022), only the drag force () is considered in this study.

*n*can be obtained by:with . represents the terminal velocity of a single particle and can be calculated by the following formula:in which

*g*is the gravitational acceleration, the density of air–water mixture must be updated by Equation (3) in every time step, and can be estimated by:

### Partial-slip boundary condition for the sediment phase

*et al.*2012; Yang

*et al.*2016). Based on Johnson and Jackson's model, Schneiderbauer

*et al.*(2012) developed an enhanced model that includes sliding and sticking collisions in one expression and provides a valuable BC for simulating the collision properties of the solid phase at the walls. A variant of Schneiderbauer's model is used in this paper to simulate the behavior of sediment particles at the wall boundary:where , , and represent the particle-wall slip velocity, wall friction coefficient at the wall, and particle-wall restitution coefficient, respectively; is the sediment phase temperature; with and ; is the normalized slip velocity and is tangential restitution coefficient. The function can be estimated as follows:

*et al.*2012) is adopted:

### Numerical scheme and procedure

For the simulations performed in this study, the Euler implicit scheme is selected for time discretization. For spatial terms: Gauss vanLeer for ; Gauss linear for other divergence terms; Gauss linear for gradient terms and Laplacian terms. The interpolation scheme is set to linear interpolation. The PCG solver and the smooth solver are chosen to solve the equations of the pressure field and the velocity field, respectively. The tolerance of each iteration of every variable is set to be less than 1.0 × 10^{−6}. Except for the sediment phase velocity, other velocity boundary conditions are set as no-slip. The pressure BC of walls is set as a fixed flux pressure boundary. The boundary conditions of walls for other scalars (, and ) are all set to zero gradient.

The governing equations are embedded into the PIMPLE scheme to form a new multiphase flow solver. The new overall solution algorithm for the proposed multiphase flow model is summarized as follows:

(Step 1) Determine and input the constant parameters in the simulation, such as gravity vector, sediment particle size, sediment density, tension coefficient, etc.

(Step 2) Initialize the fields and set the boundary conditions.

(Step 3) Use the CLSVOF solver for the new and fields by Equations (1) and (2).

(Step 4) Update the physical properties of the fluid (including air and water) phase in the whole computational domain by Equation (3) and calculate the surface tension term.

(Step 5) Compute and update sediment intergranular pressure , sediment kinematic viscosity , drag force , viscous stress , and (

*k**=**f, s*).(Step 6) Solve Equation (5) for the intermediate predicted fluid (including air and water) velocity.

(Step 7) Re-solve Equation (5) for the fluid pressure and then correct the final fluid velocity based on the new pressure field solved in this step.

(Step 8) Solve Equations (6) and (22) for the sediment phase velocity and correct the velocity at wall boundaries.

(Step 9) Solve Equation (4) to obtain the new sediment concentration distribution.

(Step 10) If < 0, set = 0; if > , set = 0.995.

(Step 11) Repeat Step 3 to Step 10 until the end of the simulation.

## MODEL VERIFICATION

The simulation of a three-dimensional silted-up dam-break flow is employed to demonstrate the solver's capability. In the experiment carried out by Duarte *et al.* (2011), the initial sediment column (1.5 m in length × 0.42 m in width × *h _{s0}* in height) was set in a reservoir (1.5 m in length × 0.42 m in width ×

*h*in height). In their study, eight sets of experiments were conducted in a rectangular tank with the size of 5.5 m in length × 0.42 m in width × 0.42 m in height. Two sets of their experiments are taken in this study for simulation: (1)

_{w0}*h*= 0.41 m,

_{w0}*h*= 0.43 m, = 1,565 kg/m

_{s0}^{3}, and (2)

*h*= 0.40 m,

_{w0}*h*= 0.22 m, = 1,671 kg/m

_{s0}^{3}. Since the initial volume fraction of sediment was not specified in their experiments (Duarte

*et al.*2011), we take the initial volume fraction as = 0.533 according to a similar experiment by Wang

*et al.*(2017). The density and viscosity of water are 1,000 kg/m

^{3}and 1.0 × 10

^{−6}m

^{2}/s, respectively. The density and viscosity of the air phase are 1 kg/m

^{3}and 1.0 × 10

^{−5}m

^{2}/s, respectively. The model parameters and material constant values are summarized in Table 1. Among them,

*K*= 1.0 × 10

^{5}Pa is used because an extremely high Young's modulus causes numerical instability and waste of computational cost (Lee 2021; Meng

*et al.*2022). The maximum Courant–Friedrichs–Lewy (CFL) number in high-concentration regions should be less than 0.005, because a too high CFL number will lead to non-convergence results, especially in sediment siltation areas. Therefore, the time step ranges from 1.0 × 10

^{−6}to 1.0 × 10

^{−4}s in this study.

Parameter . | . | . | . | b
. | a
. | . | . | . | . |
---|---|---|---|---|---|---|---|---|---|

Value | 0.1 | 5 | 0.6 | 2 | 0.1 | 0.585 | 1.5 | 0.3 | 0.8 |

Parameter | K (Pa) | (K) | |||||||

Value | 0.48 | 0.62 | 0.85 | 0.01 | 1.000 | 3.6 | 1.0 × 10^{5} | 0.33 | 288 |

Parameter . | . | . | . | b
. | a
. | . | . | . | . |
---|---|---|---|---|---|---|---|---|---|

Value | 0.1 | 5 | 0.6 | 2 | 0.1 | 0.585 | 1.5 | 0.3 | 0.8 |

Parameter | K (Pa) | (K) | |||||||

Value | 0.48 | 0.62 | 0.85 | 0.01 | 1.000 | 3.6 | 1.0 × 10^{5} | 0.33 | 288 |

*et al.*2011) with the simulated profiles of the air–water–sediment interfaces. The gate is not taken into account in this simulation, which might be the primary reason for the initial disparities between the simulated and experimental results. The comparison in the propagation of fronts (both water wave and sediment) between simulated and experimental results is shown in Figure 1(b), which shows a good agreement. It is demonstrated that the proposed Eulerian–Eulerian model is able to accurately predict the evolution of air–water–sediment interface in silted-up dam-break flows.

## NUMERICAL RESULTS AND DISCUSSION

### Geometry and conditions

*h*

_{s}represents the height of the silted-up sediment. The excellent qualitative agreement between experimental and simulated results in Section 3 implies that the simulation using the above parameters (in Table 1) is close to the real problem, except that we used finer sediment with a size of 0.0003 m and a density of 2,500 kg/m

^{3}according to the experiment by Wang

*et al.*(2017).

No. . | h_{s} (m)
. | h_{s}/H (%)
. |
---|---|---|

B1 | 0 | 0 |

B2 | 0.02 | 13.33 |

B3 | 0.03 | 20.0 |

B4 | 0.045 | 30.0 |

B5 | 0.06 | 40.0 |

No. . | h_{s} (m)
. | h_{s}/H (%)
. |
---|---|---|

B1 | 0 | 0 |

B2 | 0.02 | 13.33 |

B3 | 0.03 | 20.0 |

B4 | 0.045 | 30.0 |

B5 | 0.06 | 40.0 |

To better investigate water level fluctuations, four virtual depth probes are positioned along the center line of the channel (the *x* direction). Positions of these virtual probes are shown in Figure 2 and Table 3. In addition, four pressure probes are set to analyze the dynamics of the silted-up dam-break problem. Figure 2 (right) and Table 3 show the probe layout. To monitor the dynamic pressure during sediment collapse, the first centerline probe (P1) is positioned directly at the bottom of the channel. The other three pressure probes of the centerline are arranged on the downstream wall and the corresponding locations are 10 mm (P2), 20 mm (P3), and 40 mm (P4) above the channel bed.

Probe . | Position (m) . | Probe . | Position (m) . |
---|---|---|---|

H1 | (0.1,0.05,0) | P1 | (0.1,0.05,0) |

H2 | (0.2,0.05,0) | P2 | (1,0.05,0.01) |

H3 | (0.4,0.05,0) | P3 | (1,0.05,0.02) |

H4 | (0.6,0.05,0) | P4 | (1,0.05,00.4) |

Probe . | Position (m) . | Probe . | Position (m) . |
---|---|---|---|

H1 | (0.1,0.05,0) | P1 | (0.1,0.05,0) |

H2 | (0.2,0.05,0) | P2 | (1,0.05,0.01) |

H3 | (0.4,0.05,0) | P3 | (1,0.05,0.02) |

H4 | (0.6,0.05,0) | P4 | (1,0.05,00.4) |

### Simulation results

The present model has successfully simulated typical dam-break structures, such as sediment being deposited upstream, suspended sediment moving downstream, and sediment being exposed to air. The sediment front interacts with the surrounding water as the sediment column slides, making the water head more complicated than the initial stage. The silted-up sediment column and water flow are sheared by gravity in the early stages, which explains why the collapse occurs at the same time. The turbulence in the water flow is strong enough to suspend the sediment particles in the middle and later stages as the flow develops, causing a small amount of separated sediment particles to go downstream. However, the accumulated sediment does not move significantly, which is quite different from the suspended sediment.

### Effect of silted-up sediment height

#### Initial collapse

#### Hydrodynamics

Similarly, the measured peaks at P2–P4 suggest that the dam-break wave has reached and impacted the downstream wall. The main difference among different cases appears at the pressure peak, including the peak response time and the peak value. In clear-water dam-break flow (B1), the peak pressure value is always the highest. In silted-up cases, as the silted-up sediment height increases, the peak value of the pulsating pressure also increases.

#### Wave propagation and sediment transport

In addition, as shown in Figure 7(b), in all the cases, the propagation of sediment front over time is also traced while the reverse wave is neglected. The trajectories are relatively steep before it reaches the downstream wall at about 0.6 s, and then gradually become gentler beyond that time. As the silted-up sediment height rises, the slope of the sediment propulsion curve becomes milder. Furthermore, the time for the sediment to reach the downstream wall also differs. On the whole, the sediment propagation curves of B3 and B4 are relatively similar, but they are observably different from those of B2 and B5. It is indicated that the silted-up sediment height has a significant influence on its propagation. Theoretically, different silted-up sediment heights determine whether the dam-break process is dominated by water–sediment drag force or by stress among sediment particles. This is the decisive factor that affects the speed of dam-break flow and wave propagation. Only the propagation of the positive wave is considered in this study. The reflection and attenuation of positive and reverse waves is an interesting but quite challenging topic, which deserves our attention in the future.

Figure 8 also shows the final deposition state (the results between about 1.4 and 2.4 s) of upstream silted-up sediment for cases B2–B5 if the downstream wall boundary is ignored. This is because, after about 1.4 s, the longitudinal features of the sediment deposited upstream no longer change until the reverse wave reaches the upstream at about 2.5 s. When the reverse wave reaches the upstream, the reverse wave washes the deposited sediment, thereby changing its deposition feature. It can be concluded that in silted-up dam-break flows, the silted-up sediment height will have a significant impact on the sediment collapse process and final deposition feature.

## CONCLUSION

Based on the Eulerian–Eulerian multiphase model, a three-dimensional numerical study has been conducted to investigate silted-up dam-break flows with various silted-up sediment heights and ambient water depths. The improved CLSVOF approach and the proposed two-liquid method are used to calculate both the air–water interface evolution and sediment transport. In addition, by coupling the dynamic mechanism of sediment particles with CFD, rheology-based constitutive equations are used to simulate the deposition and erosion of sediment. The results of benchmark cases of submarine sediment collapse demonstrate that the proposed model can effectively simulate evolutions of both air–water and sediment interfaces, suggesting that the solver can reproduce gravity-dominated motion in multi-interface capturing problems. Then, the effects of silted-up sediment height on wave propagation, dynamic pressure loads, sediment transport, and sediment deposition are numerically investigated by the proposed multiphase flow model. The focus of the analysis is on the effects of silted-up sediment height on wave propagation, dynamic pressure loads, sediment transport, and sediment deposition. As the height of silted-up sediment height grows, more gravitational potential energy is consumed by intergranular stresses of the sediment phase, thereby taking longer for the dam-break flow to reach the downstream wall. However, the reverse wave reaches the upstream wall earlier. In terms of dynamic pressure loads on the downstream wall, the peak pressure values in clear-water dam-break flows are always larger than those in silted-up cases. However, in the silted-up cases, the peak value of the pulsating pressure rises with increasing silted-up sediment height. For the silted-up sediment upstream, as the silted-up sediment height declines, the time when the leftmost sediment begins to drop that indicates the beginning of a full collapse will be postponed. It demonstrates that raising the silted-up sediment height accelerates the collapse. The silted-up sediment height also has an influence on the sediment collapse process and final deposition appearance. This study provides a new insight into the emerging silted-up dam-break flow problem by proposing a new numerical model. However, the model still has some shortcomings, such as high computational time cost and less stability than a two-phase flow solver. The parameters of the model are more complex, which need to be optimized. In the future, we can study the differences in using various models (such as SPH models) for simulating silted-up dam-break flows. Turbulence might also be considered in three-dimensional silted-up dam-break flows in the future, which is a challenging but interesting problem.

## ACKNOWLEDGEMENTS

This work was supported by the Natural Science Foundation of Zhejiang Province, China (No. LY20D010009, No. LHY22E080004), Open Fund of State Key Lab of Hydraulics and Mountain River Engineering, China (No. SKHL2015), Open Fund of Key Laboratory of Flood & Drought Disaster Defense, the Ministry of Water Resources, China (No. 2021-01), National Natural Science Foundation of China (No. 51779162), National Natural Science Foundation of China (No. 51979178), National Natural Science Foundation of China (No. 52122904).

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.