Research on water flow characteristics and stress around bridges during floods

Climate change often leads to extreme rainfall, greatly increasing the probability of excessive floods in rivers (Merz et al. 2021). Many important structures in rivers, especially aged or poorly designed bridges, embankments, and other buildings, exhibit vulnerability in the face of extreme floods, often leading to safety issues for life and property. Therefore, it is imperative to have a deeper understanding of their bridge dynamic response under flood conditions and the complex water flow patterns around them. Studying the interaction between river structures and flood flow, as well as the complex flow patterns around them, is of great significance for accurately assessing the stress on structures. The damage of floods to buildings mainly stems from their strong impact force. When a significant water level difference generates high water pressure, these forces are converted into strong water flow, which directly impacts the structure and causes damage or tilting deformation. Given the severity of flood-induced structural damage, scholars worldwide have devoted considerable research efforts to uncover the underlying scientific principles. By constructing mathematical models, conducting physical model experiments, and leveraging advanced numerical simulation techniques, they have continually explored the mechanisms of flood–structure interactions and the intricate flow field variations around structures. These studies have not only enhanced our understanding of flood damage mechanisms but also provided invaluable theoretical foundations and technical support for improving the flood resistance of structures.

Many scholars at home and abroad have conducted relevant research on the structural changes caused by excessive water flow in rivers (Kyriakopoulos et al. 2022; Abdulkathum et al. 2023). Ling et al. studied the evolution of dam-break floods and their impact on structural pressure. It was found that when the flood first contacts the structure, there is a significant instantaneous impact pressure, which may reach 1.5–3.0 times the maximum pressure after the initial impact. This instantaneous impact pressure is directly proportional to the initial water depth of the reservoir (Peng et al. 2021). In this study, three-dimensional numerical simulations were performed for three typical laboratory experiments, focusing on hydrodynamic effects. Two turbulence closure models based on Reynolds-averaged Navier–Stokes (RANS) k–ε and large eddy simulation were selected. The volume of fluid (VOF) method was employed to track the free surface. The numerical model provided satisfactory results for water depth and hydrodynamic load variations over time, which were in good agreement with experimental data. The study revealed that the time-varying behavior of impact forces is significantly influenced by the upstream initial depth, with the maximum impact force exhibiting a quadratic relationship with the initial depth (Peng et al. 2023). Postacchini et al. conducted experimental tests by replicating a 1:10 scale masonry compartment, and the effect of water flow impacting the structure was obtained by moving the structure through still water. Pressure sensors were used to record the fluid-induced pressures on the four walls of the structure. Their findings indicated that the overpressure acting on the structure depends on the different flow characteristics of the front, side, and rear walls (Anfuso et al. 2021; Brocchini et al. 2022). Muñoz et al. proposed a three-dimensional non-hydrostatic RANS model using the VOF method to simulate dam-break flows. They introduced a method for recalibrating two-dimensional models to improve consistency with three-dimensional model predictions. Simulation results also showed that strong three-dimensional effects occur in regions with high river curvature, near sudden contractions, and obstacles (Munoz & Constantinescu 2020). Hien's research demonstrated that the numerical values of impact forces on obstacles from dam-break flows are more realistic in three-dimensional solutions from VOF simulations compared to two-dimensional solutions from shallow water equations. This highlights the importance of considering three-dimensional effects in simulating flood flows (Hien & Van Chien 2021). Nan et al. developed a coupled method based on the VOF method and structural dynamics to gain insights into the structural dynamic responses of bridges. Their results showed that the water flow velocity near the bridge rapidly decreased from 6 to 1 m/s, leading to bridge damage and structural deformation (Nan et al. 2023).

In summary, significant progress has been made in the numerical study of flood–structure interaction through model experiments and commercial software, providing guidance for designing structures to resist floods. However, in practical engineering cases and full-scale flood scenarios, especially under the complex hydraulic characteristics around bridge piers, there are still shortcomings in the in-depth exploration of bridge dynamic response and flow field dynamics. Recognizing this limitation, the present study builds upon the foundations of previous research by innovatively employing the OpenFOAM platform integrated with the isoAdvector method for capturing two-phase flow interfaces (Díaz-Ojeda et al. 2019; Darwish & Moukalled 2021). This approach aims to achieve precise simulation of water surface variations during flow processes, a technological breakthrough that enables more accurate capture of the nuanced details of water–structure interactions during floods. Furthermore, to enhance computational accuracy and stability, this study utilizes the snappyHexMesh tool to generate high-precision structural meshes, ensuring high resolution and reliability in the numerical simulation process. This method not only optimizes computational efficiency but also significantly improves the fidelity of simulation results to actual physical phenomena. Focusing on river velocity variations and hydraulic evolution processes under flood conditions, this research conducts detailed numerical simulations and analyses. The ultimate goal is to provide scientific evidence and practical guidance for the design, reinforcement, and optimization of river structures (such as bridges) in flood prevention efforts.

This section delves into the specifics of the numerical methods employed in our flood research. The open-source, two-phase incompressible flow solver ‘interFoam’ is adopted to simulate the flow characteristics during flood events within river channels. All descriptions in this section are based on the interFoam solver and the modifications made to enhance interface capturing. The interFoam solver utilizes the finite volume discretization method to solve the governing equations of two-phase flow (Schmitt et al. 2020).

The finite volume method, implemented in interFoam, divides the computational domain into a series of control volumes, each representing a finite volume of fluid. Conservation laws, such as the mass and momentum conservation principles, are then applied to each control volume to derive the discrete equations. This approach ensures that the numerical solution accurately reflects the physical behavior of the fluid system, especially crucial in capturing the complex interactions between water and river structures during flood events. Moreover, to improve the accuracy of interface capturing between the water and air phases, specific modifications and enhancements have been integrated into the interFoam solver. These advancements, including the adoption of the isoAdvector method for tracking two-phase flow interfaces, enable a more precise simulation of subtle variations in the water surface and the intricate hydraulic characteristics around bridge piers and other river structures.

In the interFoam solver, the spatial domain is discretized into a series of control volumes (i.e., grid cells), where the physical quantities (such as velocity, pressure, and volume fractions) are represented by their averaged values within each cell. By integrating the governing flow equations over these control volumes, a set of discrete algebraic equations is derived. These discrete equations are then solved iteratively using solution algorithms such as SIMPLE or PISO (Marsooli & Wu 2014). The SIMPLE and PISO algorithms are designed to handle the coupling between pressure and velocity fields in fluid dynamics problems. They iteratively update the pressure and velocity fields to satisfy the continuity equation and the momentum equations simultaneously. In the context of interFoam, these algorithms enable the solver to accurately capture the intricate fluid dynamics, including the interactions between the water and air phases, the effects of gravity, and any additional source terms. By discretizing the spatial domain into control volumes and applying the integral form of the flow equations, interFoam transforms the continuous fluid dynamics problem into a tractable set of discrete equations. The subsequent iterative solution process, facilitated by algorithms like SIMPLE or PISO, ensures that the numerical solution converges to a physically meaningful and accurate representation of the two-phase flow phenomenon.

This phase volume fraction α serves as a crucial indicator for tracking the interface between the two phases. In the VOF method, the continuity and momentum equations are solved over the entire domain, taking into account the effects of varying α values within each cell. This allows for the accurate simulation of the dynamic behavior of the fluid interface, including its deformation, breakup, and coalescence, which are essential features in two-phase flow phenomena.

Given the density and viscosity appearing in the continuity and momentum equations, these properties can be expressed as weighted averages based on the phase volume fraction α mentioned previously, resulting in Equation (3). The weighted averages account for the varying presence of the gas and liquid phases within each computational cell.

ρ_w represents the density of the liquid phase, ρ_a represents the density of the gas phase, and μ_w and μ_a denote their respective viscosities. These weighted averages ensure that the governing equations accurately capture the influence of the varying phase distributions within the computational domain. By incorporating the phase volume fraction α into the calculations of density and viscosity, the VOF method enables the simulation of complex multiphase flow systems with realistic phase interactions and dynamic behaviors.

Among them, c is the manually set compression factor. When c = 0, there is no compression effect. The larger the value of c, the more pronounced the compression effect.

The traditional algebraic reconstruction method simplifies the processing of ΔV_f(t,Δt) to α_fϕΔt, where ϕ represents the grid surface flux. Subsequently, the volume fraction α_f is estimated at the center of the grid surface using certain methods, which can be considered relatively rough. In contrast, the equal advection method adopts a more refined and precise processing flow. This improved method not only explains the complex fluid dynamics on the grid interface but also strives to capture the subtle differences in volume fraction evolution with higher accuracy.

• 1. Assuming that u_f remains unchanged from t to t + Δt, this is consistent with algebraic reconstruction.

• 2. Processing α_fS_f as the area occupied by phase 1, represented by A_f(t,Δt), where A_f is a function of time.

• 3. Finding the expression for A_f(t,Δt) and integrating it over t to t+ Δt. Processing ΔV_f(t,Δt) as Equation (8) in isoAdvector.

  

  (8)

isoAdvector has processed the integration of ΔV_f(t,Δt) in more detail, resulting in higher accuracy. The main steps of using isoAdvector to solve the VOF phase equation can be described as follows.

The main steps of the VOF phase equation can be expressed as follows.

• 1. Initializing ΔV_f(t,Δt) = U_fα_f,upwindS_fΔt can meet the calculation requirements for grids that are completely within phase 1 or phase 2. However, it needs to be recalculated for grids at the junction of the two phases.

• 2. By limiting the Coulomb number Co of the interface, the interface can move no more than one grid at a time (0.5 grids are more suitable), and only the flow on the surface of the grid on the interface needs to be calculated to update the next time step. This part is actually necessary and needs to be carefully combined with the isoobserver algorithm to understand that such limitations actually reduce the discussion of many additional situations (Marsooli & Wu 2014; Wu et al. 2014; Bagherzadeh Azar & Sari 2023).

• 3. For the mesh on the interface surface (0 < α < 1), reconstructing the interface (more accurately, the isosurface) and calculating the ΔV_f(t,Δt) of the interface mesh.

• 4. Updating α_p (t + Δt).

• 5. Correcting the cases where α_p (t + Δt) < 0 and α_p (t + Δt) > 1.

In computational fluid dynamics (CFD) calculations, the quality of the mesh had a significant impact on the computational results. To address this, we employed snappyHexMesh to generate high-precision structured meshes. Within the framework of these high-precision meshes, we utilized a feature-boundary adaptation method to refine the meshes specifically at the boundaries. Initially, a structured mesh with a grid size of 3 m was created for the background grid. However, to capture more detailed flow characteristics, we refined the mesh to a grid size of 1.5 m, as illustrated in Figure 5(b). Local refinement was then applied at the interfaces between the bridge and the flow field to achieve even finer structured meshes in these critical areas. To further assess the impact of mesh density on computational accuracy, we conducted an additional refinement of the background grid from 1.5 to 1 m. Upon comparing the flow information at specific monitoring points within the flow field before and after this mesh densification, we found that the changes were minimal, with errors remaining within 2% for all monitored parameters. This result indicated that the 1.5-m background grid met the criteria for mesh independence, meaning that further increases in mesh density to 1 m would not significantly alter the computational results. Therefore, we concluded that the 1.5-m background grid was sufficient to ensure computational accuracy while maintaining efficient resource utilization.

The computational domain measures 240 m in length and 210 m in width. According to field measurements, the bridge width is approximately 210 m, ensuring that the dimensions of the numerical simulation's computational domain are consistent with the actual measurements. When the size of the computational domain is determined, selecting the fluid properties becomes vital. Based on previous research (Wu et al. 2014; Wang et al. 2020), we have selected the flow parameters for the gas–liquid two-phase system within the VOF model. The parameters for the gas–liquid two-phase system are detailed in Table 1. 

After establishing the fluid properties, it is necessary to set the water velocity at the boundaries of the computational domain model. According to field measurements, the maximum river flow velocity under normal conditions is 15 m/s, whereas, during floods, the internal river velocity surges. To allow for natural development and achieve an internal velocity of 15 m/s, we set the boundary velocity at 2 m/s. The relatively shallow water depth within the river channel can be simplified for the flood process, with the water depth gradually increasing over time as the flood progresses. This entire process – establishing the computational domain, setting fluid properties, and configuring boundary conditions – constitutes the CFD approach to solving river flow problems. Assuming the grid meets the necessary requirements, the simulation is carried out using the interFoam solver within the OpenFOAM code.

Figure 7 illustrates the significant and progressively widening impact of the bridge piers on the adjacent water flow over time. The underlying cause of this phenomenon is the ‘velocity retardation effect’ triggered by the presence of bridge piers, which causes the water flow to slow down as it passes by. In accordance with Bernoulli's principle, when the flow velocity decreases, the water pressure in the vicinity correspondingly increases, resulting in the water body spreading outward. The most evident consequence is the noticeable elevation of the water surface both upstream and downstream of the bridge piers, creating the so-called backwater effect. The streamline alterations shown in Figure 8 further support this observation, with the flow velocity significantly decreasing in proximity to the bridge piers and the streamline patterns exhibiting vortex-like distortions.

The trends presented in the data in Figure 11(b) clearly demonstrate that the presence of bridge piers has a significant impact on the flow and water levels of the surrounding water bodies, creating a unique hydrological environment. Firstly, the water level at the location of the bridge piers gradually rises from a low of −20 m to stabilize at −15 m, marking a significant change of 5 m. This reflects the phenomenon of local backwater at the water level caused by the obstruction of the flow due to bridge piers. This backwater not only highlights the direct blocking effect of the bridge piers on the flow but also reveals the intricate flow patterns and vortex structures generated around the piers, which further influence the distribution of water levels.

At the intermediate interface between two bridge piers, the water level rise is even more pronounced, lifting from the same starting point of −20 to −16.9 m, which is 3.1 m higher than the water level directly influenced by the piers. This difference indicates that within the narrow channel formed by the two bridge piers, the flow velocity increases. According to Bernoulli's principle, the increase in velocity leads to a decrease in static water pressure, creating a stronger pressure difference between the banks and the space between the piers. This pushes the water body toward the center of the channel, further elevating the water level.

Additionally, Figure 12 depicts the variation in flow velocity across a mid-section of the river channel. By comparing the flow velocity data at different locations and times, it becomes evident that the flow velocity significantly decreases both upstream and downstream of the bridge piers, while the areas farther away from the piers gradually restore to normal flow velocities. This notable difference in flow velocity further corroborates the retarding effect of the bridge piers on water flow speeds and the resulting phenomenon of water level elevation.

Through a detailed analysis of the vorticity and TKE variation curves at different time instances in Figure 16, we can clearly observe the growth trends of these two physical quantities within the entire river channel over time. Notably, at the bridge piers, both vorticity and TKE reach their peak values, highlighting the significant impact of the piers on the hydrodynamic characteristics of the flow. In summary, when floods surge, bridges and their surrounding areas exhibit pronounced dynamic responses. Specifically, as the water flow encounters obstructions near the bridge piers, its velocity sharply decreases, leading to a dramatic increase in the impact force borne by the piers. This extreme force poses a high risk of damage or destruction to the bridge structure. Throughout the flood flow process, as the flood impacts the bridge, the force experienced by the bridge undergoes intense fluctuations and ultimately reaches a peak state, where the bridge's stability is at its lowest and the risk of damage significantly increases. Notably, even after the peak state, the force on the bridge remains at a relatively high level, further exacerbating the potential for damage. Concurrently, as the water flow approaches and bypasses the bridge piers, not only does its velocity undergo significant changes, but it also triggers a sharp increase in vorticity and TKE. These vortices and turbulent energy, after reaching their maximum values at the piers, gradually diffuse into the surrounding water areas, affecting a broader region. The above research results all indicate that bridges significantly affect water flow and exacerbate turbulence.

On this basis, further discussion is conducted on how to reduce the interaction between bridges and water flow, thereby maintaining the stability of bridges during flood processes. To illustrate how to reduce the interaction between water flow and the bridge pier, especially in terms of the forces acting on the bridge pier and the influence on the flow velocity, we will set different contact angles between the water flow and the bridge pier, as shown in Figures 17 and 18.

This study simulates the interaction between water flow and bridge structures during flood processes, as well as their impact on bridge performance. By monitoring key parameters such as flow velocity changes, water surface elevation fluctuations, vorticity enhancement, and turbulent energy evolution, the complex hydraulic characteristics and structural stress states of the river channel during flood processes were comprehensively revealed, as well as the influence of the interaction between water flow and bridge pier contact surfaces. Based on the above analysis, the following primary conclusions are drawn:

• (1) Structures, particularly bridge piers, are subjected to significant dynamic impacts during riverine floods. As floodwaters flow past these structures, their stress states fluctuate continuously and reach peak forces exceeding 6,000,000 N at specific moments. This sustained high-intensity force poses a severe threat to the stability of structures, highly likely leading to damage or other adverse situations.

• (2) When floodwaters traverse structures, the flow velocity undergoes marked attenuation, whereas vorticity and TKE exhibit significant increases. Over time, the influence range gradually expands, indicating that the obstruction effect of bridge piers on flood flows is not only prominent but also persistent. The enhancement of vortices and turbulence further intensifies the scouring and erosion of structures by water flows, elevating the risk of structural failure.

• (3) By adjusting the contact between the water flow and the bridge pier, the force exerted by the water flow on the bridge can be reduced, and the impact of the bridge on the water flow can be minimized.

Thank you to Professors Xiaoyu Liu and Chao Yuan for their assistance during the paper revision period. We would like to express our gratitude.

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

The authors declare there is no conflict.