Numerical simulation of submerged flow bridge scour under dam-break flow using multi-phase SPH method

1.   Introduction

Submerged flow bridge scour occurs when a bridge superstructure elevation is insufficiently high so that the bridge superstructure becomes a barrier to the flow, causing contact of the flow with the superstructure. This flow only occurs in extreme flood events or dam-break flows when the water flow is so high that it begins to approach the elevation of the bridge deck. A bridge is usually designed with the assumption of an open channel flow condition, but the water flow regime can switch to a pressure flow if the downstream edge of the bridge superstructure is partially or totally submerged during a large flood or under a dam-break flow. It is required to evaluate the effect of the submerged flow bridge scour under a bridge superstructure. In Vietnam, climate change frequently results in floods which may cause damages of bridges (Figure 1).

Investigations on submerged flow bridge scour have been reported in the studies of Abed ^[1], Jones et al. ^[2], Arneson ^[3], Umbrell et al. ^[4], Lyn ^[5], Guo et al. ^[6], Hahn and Lyn et al. ^[7], Zhai ^[8], and Kumcu ^[9]. Most of these studies deal with the experimental modeling based on Buckingham theorem. Despite significant advances in the experimental modeling, they have their own limitations in having high cost and constraints of scaling effect in the range of applicable parameters. These limitations make the experimental modeling more suitable to estimate the balance scour depth. On the other hand, the deficiency with the experimental modeling, along with the rapid advancement in computation capacity, led to the popularity of the numerical methods in the simulation of local scouring.

In the application of the numerical methods, two disciplines usually have been employed to capture the scour profiles, and they are so-called single-phase and two-phase models. The single-phase models are usually not able to consider the interparticle interactions in modeling of the scour process beneath submerged structures explicitly. However, the presence of two phases, the combination of interfacial and free-surface flows, in addition to particle entrainment of the sediment by the fluid phase are the main challenges of traditional mesh-based methods such as Finite Differences, Volumes and Elements Methods (FDM, FVM and FEM, respectively) to simulate erosion and sedimentation processes. Some inherent difficulties in some aspects of mesh-based methods can limit their application to problems that involve highly non-linear deformation such as free-surfaces and fragmentation and interior domain boundaries similar to interfacial flows ^[10]. The aforementioned limitations can be observed with multi-phase free-surface flows where the deformation of the interface is non-linear and usually fragmentation occurs in violent hydrodynamic flows ^[11].

The SPH method was developed in 1977 by Monaghan and Ringgold ^[12] and Lucy ^[13] to solve problems of astrophysics. In recent years it has been applied to other domains such as impacts of solids ^[14,15,16], fluid mechanics ^[17,18,19] and so forth. The state of the art particle methods for coastal engineering are also discussed in Ref. ^[20]. SPH has been extensively studied and extended to simulate complex problems such as multiphase solid-liquid flows ^{[21,22,23,24,25,26]}, multi-physics on interactions between fluid flow and porous ^{[27,28,29,30]}, fluid flow interactions with rigid structures or deformable structures ^[31,32,33], and coastal/ocean engineering problems. With the Lagrangian meshfree features, the SPH method has become an excellent candidate to solve several important coastal/ocean engineering problems that are not easy to be precisely/flexibly targeted by conventional related coastal/ocean engineering ^{[34,35,36,37,38,39,40]} such as simulation of floating bodies, coastal hydrodynamics, and coastal sediment transport. The SPH method continues to be developed in terms of ensuring accuracy, stability, and computation efficiency based on two main methods of two categories of weakly compressible and incompressible ones. A number of new studies have been developed to provide spatial connectivity and continuity in between calculated sources at the moving particle positions in SPH simulations such as the Background Mesh (BM) scheme ^[41] and the Moving Particle Semi-implicit (MPS) ^[40,42,43]. On the other hand, coupling of smoothed particle hydrodynamics and discrete element method (SPH-DEM) has shown the potential to be an effective numerical method. The SPH is coupled with the DEM method to simulate: 1) the formation process of the ripples resulting from a sediment transport around a swash seabed and the ripple formation from tracking individual sand particles ^[44], 2) the effect of pore water on the formation process of rolling grain ripples and in fluidization of sediment particles ^[40,45], and 3) an agitated tubular reactor and a rotating drum, showing its capability in handling complex engineering problems involving both free-surfaces and particle-fluid interactions ^[46]. An incompressible SPH method combined with the Large Eddy Simulation (LES) model can be a capable and effective tool for analyzing turbulence and eddy vortices during wave-breakwater interactions ^[47,48].

The present focus of this study lies with the bridge pressure flow scour and sediment transport of granular material of non-cohesive nature under bridge super-structure agitated by dam-break flows. The surface failure of the sediment, the rheological behavior of the scour layer and resuspension of the sediment is examined from a continuous macroscopic approach that is well suited to rheological models within SPH formulation. A code was created to solve a pressure flow scour problem based on the SPH method. The paper also presents a detailed procedure for calculating the maximum scour depth, scour profile, pressure, and shear stress of the flow under a bridge deck. Additionally, the influence of parameters such as the bridge opening height and density of the sediment on the development of the scour profile with respect to time was also numerically investigated.

2.   Multi-phase liquid-sediment SPH model

2.1.   Principals

In the SPH method, a set of arbitrarily distributed particles (points) is employed to represent the state of the system and record the movement of the system. Each particle is associated with physical parameters as mass, density, velocity, pressure, and temperature. The evolution of the system can be calculated based on the interaction of these particles and external forces. The first step to form SPH Eqs is to estimate the kernel function. For a function f at a point, whose position vector is r, it can be approximated by an integral interpolant (Eq 1).

where W is a kernel function (or smoothing function) and h is the smoothing length.

The interpolated value of function f (Eq 1) at position r can be approximated by:

where the summation is over all the particles j within the region of compact support of the kernel function, m_j and ρ_j are the mass and the density of the particle j.

The kernel function W is usually chosen to have finite support to limit the number of nearest neighbor particles involved in the summation, and thus limiting smearing effects. Several expressions of the kernel function W are used in the literature as the Gaussian function that is given by Gingold and Monaghan in 1977 ^[12], and the function of Lucy ^[13]. In this study the authors use the expression of the cubic spline of Gingold and Monaghan ^[12] that is most widely adopted in the SPH literature (Eq 3):

where α_d is a normalization factor, α_d is 1/h in 1D, 15/7πh² in 2D and 3/2πh³ in 3D, r is the distance between the particle i and j, and h is the smoothing length. In the SPH method, the length h of smoothing is defined as the radius h of the sphere of influence around the particle. The Eq for evolution of h in 2D is given by Eq 4 ^[14].

2.2.   Governing equations and SPH formulation

The continuity and momentum Eqs in Lagrangian form can be written as

where the subscript i refers to the interpolating particle and j refers to neighboring particles, $v_i^\alpha$ is the velocity vector (with relative velocity given by $v_{{\rm{ij}}}^\alpha = v_i^\alpha - v_j^\alpha$), $x_i^\alpha$ is the position vector, ρ_i is the density, m_j is the mass, $g_i^\alpha$ is the gravitational acceleration and W_ij is the smoothing or kernel function. $\sigma _i^{\alpha \beta }$ is the total stress tensor in a fluidic approach that can be written as the isotropic pressure p and the viscous stresses τ.

where δ^αβ is Kronecker's delta function, which equals 1 for α = β and 0 otherwise. In the SPH method the fluid pressure is given by the Eq of state:

where c_s is the speed of sound at the reference density, ρ₀ is the reference density, γ is a parameter chosen as γ=7 ^[49], and τ^αβ is the viscous stress. In order to distinguish between the water and sediment phases, each one requires different treatments as described in the following sections.

2.3.   Fluid phase: Constitutive equation for the water

The viscous stress tensor can be calculated from the Newtonian constitutive Eq that relates the strain rates to the viscous stresses by ^[23].

where μ represents the dynamic viscosity of water and $\dot \varepsilon _i^{\alpha \beta }$ is the deviatoric strain rate, defined as

where i and j refers to neighboring particles, Greek superscripts α and β are free indexes, γ is a dummy index and d^αβ is Kronecker's delta function. The SPH approximation may be written as

2.4.   Sediment phase: Constitutive equation for the sediment

The sediment phase is considered to be fully saturated and is modelled as a slightly compressible pseudo-Newtonian fluid, the rate dependent isotropic Newtonian fluid expression for the viscous stresses is written as. The viscous term of Eq 7 is written as

The apparent viscosity of the sediment is calculated using the Herschel-Bulkley-Papanastasiou (HBP) model ^[23,24].

where the term II_D is the second invariant of the strain rate tensor defined as

where m controls the exponential growth of stress and n is a power-law index that enables simulation of shear thinning or shear thickening behavior. Note that when m = 0 and n = 1 the model reduces to a Newtonian model, whereas when m→∞ and n = 1 a simple Bingham model is recovered. The parameter ${\tau _c}$ is the yield stress that should be defined by a yield criterion ^[23,24].

Drucker-Prager yield criterion

The yield surface prediction is modelled using the Drucker-Prager (DP) model. The DP can be written in a general form as ^[25].

where J₁ is the first invariant (trace) of the stress tensor that equals to pressure for fluids in repose - in this particular case the sediment skeleton pressure - and J_2D is the second invariant of the deviatoric shear stress tensor defined as

The parameters a and κ can be determined by projecting the Drucker-Prager onto the Mohr-Coulomb yield criterion in a deviatoric plane.

where φ is the internal friction and c the cohesion of the material. Yielding requires the square root of the second invariant J_2D of the deviatoric shear stress tensor to exceed a threshold value τ_y.

Using Eqs 15 and 19 a threshold criterion may be defined as

where τ_y is the critical shear stress that should replace τ_c in Eq 13 if the Drucker-Prager criterion is used to model the yielding mechanism of the sediment.

3.   2-D validation cases

Case 1: The SPH simulation of the sediment movement of the lat-bed sand under dam-break water is executed to compare with the model and experimental results by Zech et al ^[26]. The experimental configuration chosen for this validation is sketched in Figure 2 with the following characteristic dimensions: a water layer of depth h₀ = 0.35 m in the reservoir, and a fully saturated bed of thickness h_b = 0.125 m.). Tests were performed with bed materials: uniform coarse sand, a median diameter of 1.82 mm and a density of 2.68 ^[26]. An intense flow of water and eroded sediments is initiated after the quasi-instantaneous removal of the gate. Pictures of the flow for the flat-bed case are shown in Figure 3 at time t = 0.75 s and t = 1.25 s.

Which represents the evolution of the various levels for some characteristic times, for the lat-bed of sand material of model by Zech et al ^[26] and SPH simulation. In the model of Zech et al ^[26], various grey levels represent three regions obtained by an interface recognition procedure, while continuous lines represent the interfaces numerically computed: black lines = modeled interfaces, grey levels = experimental observations, light grey = clear water, dark-grey = moving sediment layer, and mid-grey = sediment at rest. Looking at Figure 3a and Figure 3b (time t = 0.75 s and 1.25 s, respectively), it appears that some features of the movement are well modeled of Zech et al ^[26] and SPH present simulation, such as the water surface, the scouring at the dam location (x = 0) and the sediment layer thickness. To compare the three models, experimental result of Zech et al. ^[26] and SPH present simulation, it is interesting to represent the front characteristics as well (front location X_f according to the time t) (Figure 4).

Case 2: A laboratory experiment was conducted to validate the SPH model for simulating the sediment-fluid interaction. The apparatus schematic used for the experiment is shown in Figure 5. The experimental apparatus has dimensions of 2.5 m long, 0.5 m wide, and 0.35 m high. An approximately 60 mm-thick sediment layer is placed on bottom of the flume. A sluice gate with watertight joints is lowered down to the flume bottom at the center of the test reach. The upstream water is leveled at a depth of h₀ = 0.1 m above the top of the sediment bed. The sediment bed is made of uniform coarse sand with a mean diameter of 1.5 mm. The density of the sediment was assumed to be 1800 kg/m³. A sediment trap is attached at the downstream end for collecting the eroded sediment (Figures 5 and 6). The sediment-fluid flow is captured by a digital camera (CANON D500) with a resolution of 10 Megapixels. An image analysis technique is used to capture the free surface and the sediment bed profile evolution with time during the test.

The free water surface and sediment bed profile evolution simulated using the SPH model are compared with the experimental data in Figure 7 at four different times after the dam break: 0.25 s, 0.5 s, 0.75 s and 1.0 s. Figure 8 shows the time-dependent leading edges of the water surface and the sediment bed profiles from the SPH simulation and the experimental data. The comparisons in Figs. 7 and 8 indicate a good correlation between the SPH simulation and the experimental results. However, there are a few discrepancies between the bed evolution profiles from the numerical simulation and the experiment suggesting that the SPH model under-predicted the sediment bed movement. The time duration for the wave front to reach the end of the channel was approximately 2.2 s. The simulation was continued for 5 s after the wave front reached the end of the flume in order to get adequate time for collecting a significant mass of the eroded sediment. The mass of the eroded sediment per meter wide at the end of the simulation may be calculated as dp_x·dp_y·ρ_s·N_p, where dp_x and dp_y are the particle spacing in X and Y directions, respectively, ρ_s is the sediment density, and N_p is the number of sediment particles collected at the end of the simulation. Figure 9 presents the mass of the eroded sediment per meter wide computed by the SPH simulation and the experimental results at t = 7.2 s. A spatial convergence study was performed on the initial particle spacing with dp = 0.003, 0.002, 0.0015, and 0.001 m as shown in Figure 9. Compared with the experimental result, a particle spacing of 0.0015 m is reasonably adequate.

4.   Applications

This section presents a 2D multi-phase SPH simulation of a submerged flow bridge scour under a dam-break flow. In which, the implementation of the Drucker-Prager yield criterion is used. In this numerical simulation, an upstream reservoir filled with water particles was created. The initial thickness of the sediment layer is 0.5 m, whereas the thickness of the downstream water above the sediment is 0.25 m. The characteristic parameters of the sediment in this simulation are listed in Table 1.

The schematic of the computational domain is shown in Figures 10 and 11. The height and width of the upstream water is H_w = 1.5 m and B_w = 20 m, respectively. The level of the upstream water was chosen so that the bridge would be totally submerged in the flow after the dam-break.

A spatial convergence study in terms of the position of the wave front was performed using the following initial particle spacing of dp = 0.04 m, 0.02 m, 0.01 m and 0.005 m, corresponding to the total numbers of simulation particles of 30363; 122558; 487100 and 1948000, respectively. Figure 12 shows the position of the wave front at t = 0.6 s for all particle resolutions. Hence, a particle spacing of dp = 0.005 m is reasonably adequate in this simulation. The time of simulation is 20.0 s.

To simulate the submerged flow bridge scour processes beneath the bridge deck, a 2D numerical SPH code was created in Matlab. The simulation parameters and numerical conditions are listed in Table 1.

5.   Results and discussion

Figures 13-17 represent the scour profiles, field of velocity, and pressure distribution of the water under the bridge deck at t = 1.8 s, 5.0 s, 10.0 s, 15.0 s and 20.0 s with h_b = 0.5 m, L = 3.0 m, and ρ_s = 1800 kg/m³. The location where x = 0 is the point of the deck edge and the location where y = 0 is the bed before scour. The result shows that the shape of the scour hole changes during the time between t = 1.8 s and t = 5.0 s, while the sediment dune starts to move downstream. In this time duration, the flow separates from the leading edge of the bridge deck, creating a flow separation zone and forcing the flow in the bridge opening to contract and accelerate. The high velocity flow in the bridge opening initiates sediment movement and scour. After 5.0 s, the scour gap was adequately developed, then the shape of the scour hole slightly changes during the time between t = 5.0 s and t = 20.0 s. Figure 18 shows the evaluation of scour depths at different times. The maximum scour depths are y_s = 0.9 cm, 3.0 cm, 3.8 cm, 3.8 cm and 3.9 cm at t = 1.8 s, t = 5.0 s, t = 10.0 s, t = 15.0 s and t = 20.0 s, respectively. As depicted in Figure 13, the water flow pressure peaks at the front side of the bridge deck and at the location of the scour hole while that at the rear side of the bridge deck is much smaller due to the flow separation.

When the bed shear stress induced by the fluid exceeds the critical shear stress, the initially rested sediment particles begin to move. The reduction in the transport rate can be attributed to the decrease in the bottom shear stress near the movable bed under the bridge deck. Figure 19 shows the time-averaged shear stress near the movable bed under the bridge deck at t = 1.8 s, t = 5.0 s, t = 10.0 s, t = 15.0 s, and t = 20.0 s. At the beginning of scour, the bottom shear stress increases significantly due to the high velocity flow at the scour hole beneath the bridge deck while the shear stress remains almost constant at the upstream and downstream sides. This large shear stress results in rapid sediment transport at the locations beneath the bridge deck and within the scour hole. After 20.0 s, the fluid shear stress decreases dramatically corresponding to the minimal sediment transport rate, and the form of the scour hole remains almost stable.

A parametric study using SPH simulations was also conducted considering the variations in the bridge opening height and the density of the sediment. Figure 20 shows the scour depth with varied bridge opening height (h_b). The scour depth increases as the bridge opening height decreases. Figure 21 presents the influence of the density of the sediment (ρ_s = 1800, 2000, and 2200 kg/m³) on the scour profile. Again, it can be seen that the predicted scour depth decreases as the density increases.

6.   Conclusions

A multi-phase 2D-SPH model was created to predict the submerged flow bridge scour under a dam-break flow. The Drucker-Prager yield criterion was used in the SPH simulation in order to model the relevant features of the erosive process considering the mechanical and hydraulic properties of the sediment bed. Two validation cases show that the results from the SPH model created in this study compare well with both the results of Zech et al. ^[26] and the laboratory experimental results conducted by the authors. The SPH model shows the capability of predicting the flow velocity, flow pressure fields, and scour profiles beneath a bridge deck, as well as the progress of fluid-sediment interfacial shear stress under a submerged flow bridge due to a dam-break flow. The simulation results also show great influence of parameters such as bridge opening height and the density of the sediment on the evolution of the scour profile with respect to time.

Conflict of interest

The authors declare no conflicts of interest.