The RNG k–ε turbulence model is based on the improvement and modification of the standard model. The RNG k–ε adds an additional correction term to the ε equation, which improves the accuracy of calculating the rapid flow conditions.

In this study, we used the RNG turbulent flow model and the DPM discrete phase of FLUENT 16.0 software for numerical calculations. VOF was used for free interface tracking (Thinglas and Kaushal, 2008). The continuity equation, momentum equation and K, and the ε equation of the RNG turbulent flow model can be expressed as follows: ∂ ρ ∂ t + ∂ ρ u i ∂ x i = 0 , (1) ∂ ρ u i ∂ t + ∂ ( ρ u i u j ) ∂ x i = − ∂ ρ ∂ x j + ∂ ∂ x j [ ( μ + μ i ) ( ∂ u i ∂ x j + ∂ u j ∂ x i ) ] + ρ B i , (2) ∂ ( ρ k ) ∂ t + ∂ ( ρ u i k ) ∂ x j = ∂ ∂ x i [ ( μ + μ i σ k ) ∂ k ∂ x i ] + G k − ρ ε , (3) ∂ ( ρ ε ) ∂ t + ∂ ( ρ u i ε ) ∂ x j = ∂ ∂ x i [ ( μ + μ i σ ε ) ∂ ε ∂ x i ] + C 1 ε ρ ε 2 k , (4) where ρ is the average density of the volume fraction, μ is the molecular viscosity coefficient, t is the time, x_i is the spatial coordinate in the i direction, u_i is the velocity component in the i direction, B_i is the force per volume, G_k is the turbulent kinetic energy induced by the mean velocity gradient, k is the turbulent kinetic energy, ε is the turbulence dissipation rate, and σ ε and C 1 ε are the constant values used in the turbulence model ( σ ε   =   0.718 , C 1 ε   =   1.68 ).

The orbits of discrete-phase particles (Akhtar et al., 2007) were solved in FLUENT by integrating the differential equation for the forces acting on the particles in the Rasch coordinate system. The equilibrium equation for the forces acting on the particles (particle inertia = various forces acting on the particle) in the Cartesian coordinate system had the following form (x-direction): d u p d t = F D ( u - u p ) + g x ( r p − r ) r p + F x , (5) F D = 18 m r p d p 2 C D Re 24 , (6) Re = r d p | u p − u | m , (7) where F_D represents the mass traction of particles, u represents the fluid phase velocity, u_p represents the particle velocity, g_x represents the gravity acceleration in the x direction, m is the hydrodynamic viscosity, r represents the fluid density, r_p is the particle density, d_p represents the particle diameter, and Re represents the relative Reynolds number.

The traction coefficient C _D can be expressed as follows: C D = a 1 + a 2 Re + a 3 Re . (8) For spherical particles, a₁, a₂, and a₃ in the above equation are constants for a range of Reynolds numbers.

The other forces included in the force balance equations for particles may be important in some cases. The most important of these “other” forces is the “apparent mass force” (additional force), which is the additional force caused by the acceleration of the fluid around the particle, and it is expressed as follows: F x = 1 2 r r p d d t ( u − u p ) . (9) When r > r p , the apparent mass force cannot be ignored. The additional forces due to the fluid pressure gradient present in the flow field are as follows: F x = r r p u p ∂ u ∂ x . (10) According to the soil particle classification of test soils derived from the U.S. Department of Agriculture (USDA), the Loess Plateau has a relatively large proportion of silty soil (0.02–0.002 mm). The particle diameter of the simulated soil was set at 0.02 mm using homogeneous sediment particles.

The outlet section is assumed to be a fully developed turbulent flow. The inlet section is a given water-level value, and the turbulence energy K is used in the empirical formula. The turbulence intensity I is calculated using the following equation: I = 3 / 2 ( Re ) − 1 / 8 , (11) k = 3 2 ( u ¯ I ) 2 , (12) ε = C μ 3 / 4 k 3 / 2 L , (13) where u ¯ represents the average flow velocity, I is the turbulence length scale (I = 0.07 L), L represents the hydraulic diameter, and C μ represents an empirical constant ( C μ   = 0.09 ).

Sediment particles lose some kinetic energy when colliding with the wall. To account for this, we defined the geometric bottom of the model as a reflective wall. The recovery coefficient after collision indicates the kinetic energy loss of the sediment particles. We set the normal phase and tangential recovery coefficients of the particles colliding with the wall, assuming the same amount of energy is lost in each collision (Ye, 2019). The specific formulas are as follows: e n = 0.5769 − 0.00231 θ + 0.04766 θ 2 ( 0 < θ ≤ π 2 ) , (14) e t = 0.4747 + 0.4254 θ + 0.06756 θ 2 ( π 10 ≤ θ ≤ π 2 ) , (15) where e_n represents the normal phase recovery coefficient of sediment particles colliding with the wall, e_t represents the tangential phase recovery coefficient of the collision between sediment particles and the wall, and θ represents the direction of sediment particles colliding with the wall and the angle of the collision surface.

In the DPM model, the deposition rate is equal to the ratio of wall deposition mass flow rate to the area, according to the particle mass balance equation. The equation is as follows: R A = ∑ p = 1 N p a r t i c l e s m p A f a c e , (16) where R_A represents the sediment deposition rate, m_p represents the mass flow rate of the particle stream, and A_face represents the area of the wall at the particle impact boundary.

In this test, the finite volume method is used as the control equation, and the second-order implicit scheme is used in the discretization of the time item. The PISO algorithm is used to solve and control the coupling of velocity and pressure in the equation. The VOF method is used to track and simulate the free surface and two-phase flow of air and water. The free water surface is established by the geometrical reconstruction scheme.

Using drone aerial photography technology, Figure 2 is taken from the top of the stilling basin. The stilling basin is located at the Guandigou 1^# check dam (110˚37'E, 37˚58'N) in Suide County, Shaanxi Province, in the Jiuyuangou watershed. The main ditch of Guandigou is 18 km long, with an average slope gradient of 1.15%, a “V” shaped ditch profile, a gully density of 5.34 km/km², and an elevation of 820–1,180 m. The topography of the watershed is fragmented, the gullies are crisscrossed, and soil erosion is a serious issue. The average soil erosion modulus is 14,000 t/(km² a) along the middle and lower reaches of the Wuding River.

The model was built at a 1:1 scale and consisted of a steep slope, a dissipation pond, a gradual section, and a tail channel (Figure 3). As shown in Figure 3, the overflow weir connected with the steep slope was 24 m high, so the starting point of the steep slope was 24 m higher than the ground elevation, with a slope drop ratio of 1:5. The steep slope was 120 m long, the stilling basin was 11 m long and 2.2 m deep, the gradient section was 8 m long, and the tail channel was 20 m long, with a slope drop ratio of 0.025.

The model’s geometry consisted of a quadrilateral grid defined by ICEM CFD 16.0. Grid accuracy varied based on proximity to the extremes of the calculation domain. The grid near the bottom of the calculation domain was more precise (0.02 mm), while the grid near the top of the calculation domain was less precise (0.10 mm). We accounted for this variation because the steep slope had only a thin layer of water flow. The entire grid consisted of 285,100 cells (Figure 3). Its grid average mass was greater than 0.92, and no negative volume appeared. This means that the model had great mesh quality. Simulations were carried out using the Fluent 16.0 commercial package. The uncoupling arithmetic method was used to separate and solve model equations: pressure–velocity coupling was defined by the pressure-implicit with splitting of operators (PISO) algorithm. The left side of the calculation domain represented the point of traffic entry, while the right side of the calculation domain represented the pressure outlet. The top of the calculation field represented the pressure, and the bottom of the calculation domain represented the reflective wall.

The test is based on a spillway designed for a 20-year flood. The bottom width of the check dam’s spillway trapezoidal section is 10.5 m, and the side slope ratio is 1:1. The unit width flow rate of the spillway inlet is the ratio of the maximum spillway flow rate to the bottom width of the inlet section. The change in the spillway discharge flow rate is in 1.5 times increments (Table 1). To obtain the maximum flow rate of the spillway, we consulted the existing check dam design information. The calculation formula is as follows: Q D F = Q F ( 1 - V S F / W F ) , (17) where Q _DF represents the maximum flow rate of the spillway, Q _F represents the 20-year design flood flow (105 m³ s⁻¹), V_SF represents the stagnant flood volume (550,000 m³), and W_F represents the maximum flood releasing capacity of the spillway (100,000 m³).

The annual sand discharge of the reservoir does not exceed 40%. The average annual sand discharge is 20% (Xin et al., 2012; Wang et al., 2020). The erosion modulus of the Loess Plateau is 1000 t / ( k m 2 ⋅ y e a r − 1 ) . The formula for annual sediment discharge is as follows: M 0 = 0.2 0 ∗ F K = 1000 kg, (18) where M₀ is the annual sand discharge (kg), F is the check dam control watershed area (F = 5 km²), and K is the annual average erosion modulus [k = 1,000 t/(km² year⁻¹)].

Assume that the annual sand discharge mass is obtained during one rainstorm and all the sediment is discharged by the spillway. The unit width sand mass flow rate is then obtained by dividing the sand discharge mass flow rate by the bottom width of the inlet section. Different unit width sand discharge mass flow rates were defined for each working condition, including CK (clear water), CL₁ (M = 0.1 kg s⁻¹), CL₂ (M = 0.2 kg s⁻¹), and CL₃ (M = 0.3 kg s⁻¹). The median particle size of the sediment in the Yellow River field sub-high sand-bearing flood is generally 0.01–0.03 mm (Wang et al., 2020). We used homogenous sand with a grain size of 0.02 mm for the simulation.

Jiuyuangou watershed is located in a temperate semiarid region. The average precipitation is 469 mm, with rainfall mostly concentrated over 6–9 months, falling in heavy rainfall events. The main function of the check dam stilling basin is to dissipate the energy of the rising water in front of the dam during heavy rainstorms, thus protecting the downstream farmland from being washed away. Data on the check dam’s stilling basin were obtained from the stilling basin flood control center, which provided data on the sequent water depths when the flow rates were 2.5, 5, 7.5, and 10 m³/s. Because these flow rates are low, the sediment content added to the stilling basin was ignored in the calculations. Table 2 shows the measured and empirical formulas for calculating values at different flow rates, and the empirical formula is shown in Eq. 21. The error of the conjugate water depth ratio of the measured value to the calculated value is within 6%. Because the error is small, this indicates that the true value can be replaced by the calculated value of the empirical formula. Second, Figure 4A shows the simulated and empirical formula-calculated values under CK₃ processing. The horizontal coordinate is the Froude number (Fr), and the vertical coordinate is the conjugate bathymetry ratio. The errors of the simulated and empirical formula-calculated values are within 10%. This shows that the operating conditions are reliable in clear water and that these can be used for other operating conditions as well. Finally, the sediment deposition at the bottom of the pool could not be obtained accurately due to the unstable water flow. Therefore, to verify the reliability of the sediment deposition simulation, the settling velocity of the modeled ( ω m ) single-particle sediment was compared with the existing empirical ( ω c ) equation (Camenen, 2007). Camene proposed a simple, reliable, and general equation for the settling velocity of the particles that accounts for their shape and roundness. The vertical coordinate of Figure 4B shows the ratio of calculated ( ω c ) to simulated values ( ω m ), which is lower than 1.3. This indicates that the simulated results are similar to those theoretical values calculated by the empirical formula. h ′ = h ′ 2 ( 1 + 8 q 2 g h ′ ′ 3 − 1 ) , (19) h ′ = h ′ 2 ( 1 + 8 q 2 g h ′ 3 − 1 ) , (20) η = h ′ h ′ = 1 2 ( 1 + 8 F r 1 2 − 1 ) , (21) where q is the pre-jump section unit width flow ( m 2 ⋅ s − 1 ), h ' is the pre-jump water depth (m), h ' ' is the post-jump water depth (m), and η is the conjugate water depth ratio.

To calculate the energy dissipation rate (K_j) of the stilling basin at different flow rates, we compared the energy of the initial section (E₀) at the bottom of the stilling basin to the energy of the exit section (E₁) at the tail channel. The total energy of each section can be expressed as E i = Z i + p i γ + v i 2 2 g . Here, Z i is the potential energy, p i is the average pressure, and v i is the average flow velocity of the cross section. The energy dissipation rate of each section can be expressed as η = E 1 - E 0 E 0 .

The energy dissipation rate increased with increasing Fr (Table 3), and the energy dissipation rate was greater in clear than in sandy water. The reason is that the sand-bearing water flow accelerates the flow velocity of the slope surface water flow so that the flow velocity of the water flow into the stilling pool is faster. At a constant unit width flow rate, the energy dissipation rate increased as the sand discharge mass flow rate increased. This shows that under the test conditions, the sand concentration of the flow affected energy dissipation: when Fr > 4.5, K_j was greater than 50% and when Fr < 4.5, K_j was less than 50%. Figure 5 illustrates the water phase cloud at the bottom of the stilling basin for different inlet unit width flows. As the inlet unit width flow increased, the flow pattern at the bottom of the basin destabilized, thereby increasing energy dissipation. At a unit width flow rate less than 3 m²/s, the flow pattern was more stable.

Figure 6 shows the sediment deposition at the bottom of the stilling basin. The sediment was first deposited at the end of the pool bottom. Then the deposition height and length gradually increased. After reaching the maximum deposition height of 0.2 m, the deposition length continued to increase until the pool bottom was covered (Figures 6D2,D3). In addition, as the mass flow rate increased and the inlet unit width flow rate decreased, the maximum deposition rate at the bottom of the pool increased. As the inlet unit width flow rate decreased, the flow pattern of the basin bottom gradually stabilized (Figure 5). This indicates that the more stable the flow pattern of the pool bottom, the lower the Fr and the greater the sediment deposition at the bottom of the pool.

The effects of dynamic and sediment parameters on the wall deposition rate were analyzed using the Pearson correlation analysis (Figure 7). The inlet flow rate (q₀) was found to be strongly and positively correlated with the inlet velocity (v₀), turbulent kinetic energy (E), and shear force ( τ ), with a correlation coefficient as high as 0.92. The deposition rate (R_A) at the bottom of the pool was negatively correlated with the inlet flow, inlet velocity, turbulent energy, and shear, with correlation coefficients greater than 0.73. The deposition rate was positively correlated with the mass flow rate of sand discharge (M), inlet sediment concentration (C₀), and average sediment concentration (C). The correlation coefficient between sedimentation rate and inlet sediment concentration was as high as 0.97. Thus, the relationship between sedimentation rate and inlet sediment concentration was linearly fitted, with R ² = 0.94 (Figure 8). The inlet sediment concentration is the ratio of the mass flow rate of discharged sand to the inlet flow rate. The final equation for the relationship between the average deposition rate and the mass flow rate and inlet unit width flow rate is as follows: C 0 = M q 0 , y A = 69 . 75 M q 0 + 1.75 = 69 . 75 C 0 + 1 . 75 , where y_A is the deposition rate ( k g ⋅ m − 2 ), M is the mass flow rate of sediment discharge ( k g ⋅ s − 1 ), q₀ is the inlet unit width flow rate ( m 2 ⋅ s − 1 ), C₀ is the inlet sediment concentration per unit area ( k g ⋅ m − 2 ), and H₀ is the water inlet height (m).

The energy dissipation rate is an important index for measuring the effect of energy dissipation. Based on the analysis of the energy dissipation rate in a previous article, many scholars believe that Fr has a significant effect on the initial energy dissipation rate. Li et al. (2018) believes that the larger the Fr, the higher the energy dissipation rate. Zhang et al. (2017) proposes that when Fr is 4.5–9.0, the energy dissipation effect improves because the water jump is stable. Sun et al. (2019) further refined the range of the effect of Fr on the initial energy dissipation rate and concluded that when Fr is 2.55–4.5, the energy dissipation rate (Kj) ranges from 20 to 45%, and when Fr is 4.5–9, the energy dissipation rate (Kj) ranges from 45 to 85%. Here, we found that the variation in the energy dissipation rate (Kj) between Fr in the pre-leap (E0') and the post-leap section of the dissipation pool basically conformed to this law, even under different working conditions. When Fr was > 4.5, the energy dissipation rate Kj was greater than 45%. When Fr was < 4.5, the energy dissipation rate Kj was less than 47%.

Also, in Figure 9, it is shown that the higher the inlet flow, the higher the turbulent kinetic energy, the more vortices are generated, and the greater the energy dissipation rate. We further believe that the more vortices, the more significant the flow rate reduction and the better the energy dissipation effect. This is consistent with the conclusion of other experts as well. Tan et al. (2020) believed that the surface vortex zone is significantly larger than the bottom vortex zone in the second-stage force elimination pool, and that the water flow is more stable and the energy dissipation rate is larger than that in the single-stage pool. Dong et al. (2016) found that a non-complete wide tail pier force pool can produce a large transverse velocity gradient and that this transverse velocity gradient produces additional turbulent shear and lateral flow. In this way, the energy dissipation effect is improved, as desired.

The maximum lift force generated by turbulent pressure fluctuations acting on the bottom of the stilling basin results in poor sediment stability (Bowers and Tsai, 1969). If the lifting force is greater than the gravity of the sediment particles, no sediment will be deposited. If the lifting force is less than the gravity of the sediment particles, the sediment will be deposited on the bottom of the basin. The bottom surface of the water flow is gradually transformed from a relatively smooth basin to a submerged sediment topography with a certain undulating height, and the resulting additional friction at the bottom leads to a relative increase in water flow resistance (Zhu, 1982; Wei 2013).

As the inlet unit width flow rate increased, the siltation height at the bottom of the basin did not exceed 0.2 m. The effect of inlet unit width flow on the maximum silt thickness was insignificant. Chen et al. (2019) predicted the deposition of raw Yangtze River water into the basin’s transfer pond, concluding that when the shape of the inlet is fixed, the sediment deposition distribution does not change much. Guo (2017) analyzed the sediment deposition pattern of the river and determined that the maximum deposition of the river would not exceed 0.03 m, and the deposition in most areas was below 0.017 m. Our results are consistent with Chen’s conclusion, suggesting that a fixed inlet shape results in decreased significance of flow pattern on maximum sedimentation height at the bottom of the basin. Flow is the main factor affecting the sand discharge rate at the outlet of the stilling basin. The flow rate has a significant effect on the amount of sediment siltation at the bottom of the stilling basin. By simulating sediment deposition in inverted siphons, Bagchi (2012) concluded that the deposited sediment gradually decreases as the flow rate increases.