Simulation Methods for Flash Floods and Debris Flows Based on Rheological Evolution Characteristics

Abstract

This invention belongs to the field of flash flood and debris flow disaster prevention and mitigation, and provides a method for simulating flash flood and debris flow processes based on rheological evolution characteristics. The method includes the following steps: establishing a generalized numerical model of the target channel; using a CFD-DEM coupled model to simulate the transport process of sediment particles in the generalized numerical model under the action of water flow or debris flow slurry; and, based on the semi-empirical relationship between the drag coefficient and Reynolds number of sediment particles carried by flash floods and debris flows, derived in this invention, realizing real-time adjustment of rheological parameters during the simulated flash flood and debris flow transport process through the bidirectional coupling mechanism of the CFD-DEM coupled model. This invention's method can more realistically simulate the spatial distribution of sediment particles during flash flood and debris flow transport, increasing the accuracy of flash flood and debris flow process simulation, and solving the problem that existing technologies struggle to simulate the rheological evolution phenomenon of flash floods and debris flows caused by changes in volumetric sediment content.

Description

Technical Field

[0001] This invention belongs to the field of disaster prevention and mitigation of flash floods and debris flows, and relates to a method for simulating flash flood and debris flow processes based on rheological evolution characteristics. Background Technology

[0002] Debris flows, composed of large amounts of silt, rocks, and water, possess powerful erosive and destructive forces. The rheological characteristics of debris flows reflect their physical properties and kinematic characteristics, directly influencing their erosion process. Currently, research on the rheological characteristics of debris flows largely relies on experimental measurement techniques and empirical rheological models. Among existing technologies, the Bingham model and power-law model are commonly used to describe the rheological characteristics of debris flows. The key parameters of these models mainly include the yield stress and viscosity coefficient of the slurry. For example, Huang Jiangcheng et al. conducted debris flow experiments using representative soils from the Xiaojiang River basin in Yunnan and the Bailongjiang River basin in Gansu, exploring the effects of particle size distribution, clay content, and mineral composition on the rheological characteristics of the slurry. Contreras and Davies, by analyzing the rheological characteristics of debris flows carrying coarse particles, revealed the influence of coarse particles on the yield stress variation. Yu Jie et al. revealed the relationship between the fractal characteristics of convective debris flows and their kinematic parameters such as flow rate and sediment transport rate. Kaitna et al. studied the effects of sediment gradation and content on pressure distribution within debris flow slurry, explaining the role of sediment in increasing the fluidity and long-range movement of debris flows. Sun et al. used numerical simulation to explore the evolution range of debris flow processes under different pore water pressure ratios in the Izu mountainous region of Shizuoka Prefecture, Japan. Nie Yipin et al. simulated the deposition process of sediment particles of different sizes at the mouth of debris flow channels, exploring the influence of volumetric sediment content (Cv) on the distribution characteristics of sediment particles. Han et al. improved the simulation effect of debris flow erosion patterns and inundation areas based on cellular automata models.

[0003] Currently, most numerical simulation studies of debris flows, both domestically and internationally, employ fixed rheological parameters. This means that the rheological parameters of the debris flow must be preset before the simulation begins and remain unchanged during the simulation, although the parameter values ​​vary significantly between different simulations. However, the volumetric sediment load (Cv) of debris flows changes during transport, and this change in Cv leads to changes in the rheological parameters. Therefore, existing numerical simulation methods using fixed rheological parameters cannot accurately reflect the rheological evolution of debris flows caused by changes in Cv during transport. To address this, Shen et al. proposed considering the rheological evolution of debris flows caused by bed erosion, revealing the significant impact of bed erosion-induced rheological properties and topographic changes on the propagation of flowing debris flows. Meyrat et al. introduced the Voellmy mixed rheological model into their numerical model to obtain the erosion patterns and density distributions of debris flows, reproducing the conversion process between debris flows and floods or hyperconcentrated flows. However, these studies still require prior assumptions about the erosion rate of the solid phase, making it difficult to simulate the movement of eroded sediment at the front. Furthermore, while the aforementioned single-fluid models have advantages such as fewer computational parameters, lower computational cost, and ease of use, they assume that the velocities between the solid and liquid phases of the debris flow are essentially the same, neglecting the relative motion between the two phases. This assumption is clearly invalid when the sediment content of the debris flow is low. Therefore, using a discrete medium model to calculate the movement of sediment particles and coupling it with a continuous medium model to obtain the fluid motion is a feasible method for simulating debris flows.

[0004] Computational fluid dynamics-discrete element method coupled simulation (CFD-DEM) can quantify the forces between phases, capturing gully deformation and sediment movement from a particle perspective, and has been widely used in the simulation of flash floods and debris flows. Nie et al., based on the CFD-DEM method, constructed an improved Weibull function model to quantify the instantaneous sediment transport rate variation characteristics of loosely deposited sediment particles during flash floods and debris flows. Li et al., using the CFD-DEM method, found that the protective efficiency of protective nets is closely related to the Froude number of the debris flow slurry. Lei Ming et al. used the CFD-DEM method to reveal the relationship between the evolution of water level along the channel and the gully erosion and deposition deformation and sediment supply conditions in variable gradient flash floods and debris flows. Kong and Guan, based on the CFD-DEM method, studied the blocking effect of protective structures on different flash floods and debris flows. However, these studies still require pre-defining rheological parameters and keeping them constant during the simulation process. It is difficult to simulate the rheological evolution of debris flows caused by changes in volumetric sediment content (Cv), and the accuracy of the simulation needs to be improved. Summary of the Invention

[0005] To address the problem that existing technologies require predefined rheological parameters and must maintain these parameters constant during coupled simulation of flash floods and debris flows, thus failing to simulate the rheological evolution characteristics of flash floods and debris flows caused by changes in volumetric sediment content, this invention provides a flash flood and debris flow process simulation method based on rheological evolution characteristics. The method considers the impact of changes in volumetric sediment content on the rheological characteristics of flash floods and debris flows during transport. By adjusting the rheological parameters of the flash flood and debris flow in real time during the simulation, it achieves the simulation of the flash flood and debris flow process under varying volumetric sediment content conditions. This more realistically simulates the spatial distribution of sediment particles during the transport of flash floods and debris flows, increases the accuracy of flash flood and debris flow simulation, and provides more scientific guidance for the prevention and control of flash flood and debris flow disasters.

[0006] To achieve the above-mentioned objectives, the technical solution adopted by the present invention is as follows:

[0007] A method for simulating flash floods and debris flows based on rheological evolution characteristics includes the following steps:

[0008] S1. The hilly gully to be simulated for the flash flood and debris flow process is taken as the target gully. The target gully includes the interconnected debris flow formation area and debris flow flow area. The debris flow flow area gully has loose deposits. A generalized numerical model of the target gully is established based on the topographic information of the target gully and the gradation information of the sediment particles in the loose deposits.

[0009] S2, the CFD-DEM coupled model is used to simulate the transport process of sediment particles in a generalized numerical model under the action of water flow or debris flow slurry. The steps are as follows:

[0010] S21, Read the target channel generalized numerical model;

[0011] S22, the DEM module executes the DEM solving loop. In the DEM solving loop, the particle neighbor list is first read and the contact and collision between particles are identified. The force on each particle is calculated. Then the motion state of each particle is determined, the position and velocity of each particle are updated, and the position and velocity information of the particles are passed into the coupling module of the CFD-DEM coupled model.

[0012] S23, the position and velocity information of the particles obtained by the DEM module is passed to the CFD module through the coupling module of the CFD-DEM coupled model. The CFD solution loop is executed in the CFD module. In the CFD solution loop, the fluid control equation is solved and the flow field information is updated. When the flow field converges, the flow field information is passed to the coupling module of the CFD-DEM coupled model.

[0013] S24. In the coupling module of the CFD-DEM coupled model, the fluid velocity and mesh porosity at the location of each particle are obtained, the relative velocity of the particle fluid is calculated, and it is determined whether the particle is in the water flow. If the particle is in the water flow, the rheological parameters, including yield stress and viscosity, are calculated according to the porosity of each mesh. If the particle is in the air, the drag coefficient of the particle is calculated according to the DiFelice formula. For particles in the water flow, it is determined whether the particle has started. If the particle has not started, the drag coefficient of the particle is calculated according to the DiFelice formula. If the particle has started, the equivalent Reynolds number of the particle is calculated according to equation (4). Then, the separation angle and projected area coefficient are calculated through equations (11) and (13). Then, the drag coefficient of the particle is calculated through equation (17). Finally, the total drag on each particle is calculated, and the total drag information on the particle is transmitted back to the DEM module.

[0014]

[0015] In the formula, Re * ρ is the equivalent Reynolds number of the particle. f Let d be the density of the fluid containing the particle, d be the particle diameter, and v be the density of the fluid containing the particle. max η is the maximum settling velocity of particles in Bingham fluid, η0 is the viscosity, and ε * Let η be the dimensionless settling velocity, Re' be the Reynolds number of the particle, and η be the particle number. * The apparent viscosity is dimensionless, and θ is the separation angle. is the projected area coefficient of the separation zone, and Cd is the drag coefficient of the particles.

[0016] S25, repeat the operations of S22 to S24, and continue the calculation for the next time step.

[0017] In step S2 of the above technical solution, when using the CFD-DEM coupled model to simulate the transport process of sediment particles in the generalized numerical model under the action of water flow or debris flow slurry, the water flow rate or debris flow slurry flow rate is determined by analyzing historical rainfall data at the target channel.

[0018] The water flow rate can be determined by the following steps: based on the design flood control index of the area where the target channel is located, combined with historical rainfall data and hydrological analysis, determine the flood flow rate for a specific return period, use the flood flow rate for that specific return period as the original flow rate, and scale the original flow rate according to the model scale to obtain the water flow rate used in the generalized numerical model.

[0019] The debris flow slurry flow rate can be determined by the following steps: Based on historical rainfall data and hydrological analysis results of the target channel area, combined with the correlation analysis data between historical rainfall and debris flow flow rate, the debris flow flow rate triggered by rainfall at a specific return period is determined. The debris flow flow rate at that specific return period is used as the original flow rate. The original flow rate is scaled according to the scale of the generalized numerical model to obtain the debris flow slurry flow rate used in the generalized numerical model.

[0020] In the above technical solution, when step S2 uses the CFD-DEM coupled model to simulate and calculate the transport process of sediment particles in the generalized model under the action of debris flow slurry, the volume concentration of the debris flow slurry is determined by analyzing the historical rainfall data, debris flow sediment yield and scale at the target channel.

[0021] In step S22 of the above technical solution, calculating the force on each particle refers to calculating the nonlinear elastic normal and tangential forces between particles and between particles and boundaries based on Hertzian contact theory and Mindlin and Deresiewicz theory, respectively; after determining the force on each particle, the motion state of each particle is determined according to Newton's second law. This method of calculating the force on each particle and determining the motion state of each particle is an existing method.

[0022] In step S23 of the above technical solution, the fluid control equations are the fluid continuity equation and the NS equation. The PISO algorithm is used to solve the fluid control equations, and the method of using the PISO algorithm to solve the fluid control equations is an existing method.

[0023] In step S24 of the above technical solution, when determining whether a particle is located in the water flow, the density value obtained by weighted average of the fluid volume fraction of the grid where the particle is located is used to determine whether the particle is located in the water flow. This determination method is an existing method.

[0024] In step S1 of the above technical solution, a generalized numerical model of the target gully is established in the open-source software Gmsh based on the topographic information of the target gully and the gradation information of the sediment particles in the loose deposits.

[0025] Compared with the prior art, the technical solution of the present invention has the following beneficial technical effects:

[0026] 1. This invention provides a method for simulating flash flood and debris flow processes based on rheological evolution characteristics. This method is proposed based on the semi-empirical relationship between the drag coefficient and Reynolds number of sediment-carrying particles in flash floods and debris flows, derived theoretically. The method considers the influence of changes in volumetric sediment content on the rheological characteristics of flash floods and debris flows during transport. By adjusting the rheological parameters of the flash flood and debris flow in real time during the simulation, it achieves the simulation of flash flood and debris flow processes under varying volumetric sediment content. Compared to existing technologies that require pre-defining rheological parameters before coupled simulation and maintaining these parameters constant during the simulation, making it difficult to simulate the rheological evolution of flash floods and debris flows caused by changes in volumetric sediment content, the method described in this invention can more realistically simulate the spatial distribution of sediment particles during flash flood and debris flow transport, increasing the accuracy of flash flood and debris flow process simulation. This provides more scientific guidance for the prevention and control of flash flood and debris flow disasters and also contributes to in-depth research on the microscopic mechanisms of flash flood and debris flow transport in hilly areas.

[0027] 2. Based on the semi-empirical relationship between the drag coefficient and Reynolds number of sediment particles carried by debris flows during their movement, obtained through theoretical derivation, the method of this invention realizes real-time adjustment of rheological parameters during the simulated debris flow transport process through the bidirectional coupling mechanism of the CFD-DEM coupled model. This real-time adjustment helps to more accurately capture the interaction between the solid phase and the fluid phase, and provides more accurate predictions for the debris flow movement and scouring process.

[0028] 3. The method described in this invention considers the influence of sediment volume fraction on fluid rheological parameters and the changes in sediment volume fraction during the movement of flash floods and debris flows. This allows for a more realistic simulation of the spatial distribution of sediment particles during the transport of flash floods and debris flows. The method of this invention can reveal the spatial distribution law of rheological parameters of flash floods and debris flows, providing strong technical support for the prevention and control of flash flood and debris flow disasters. Attached Figure Description

[0029] Figure 1 The figure shows the relationship between the equivalent Reynolds number and the drag coefficient when particles settle in Bingham fluid. The data points in the figure represent the experimental results from reference [4], and the curves in the figure represent the theoretical calculation results of equation (17).

[0030] Figure 2 This is a schematic diagram of the experimental setup for the dam failure experiment in reference [5].

[0031] Figure 3 This is a comparison chart of the experimental data in reference [5] and the simulation results in this invention.

[0032] Figure 4 A schematic diagram of the calculation process of the method described in this invention.

[0033] Figure 5 This is a schematic diagram of a mobile bed model for flash floods and debris flows.

[0034] Figure 6 It refers to the erosion characteristics of the gully bed after contact between mountain torrents and debris flows and particles, taking into account the rheological evolution.

[0035] Figure 7 It refers to the erosion characteristics of the gully bed after the contact between mountain torrents, debris flows and mud and sand, without considering the rheological evolution.

[0036] Figure 8 It takes into account the characteristics of the movement of mountain torrents and debris flows, as well as the changes in the elevation of the gully and the water surface.

[0037] Figure 9 This refers to the characteristics of the movement of the gully bed and water surface elevation changes in flash floods and debris flows without considering rheological variations.

[0038] Figure 10 This is a simulation result of the spatial distribution of sediment volume fraction in flash floods and debris flows considering rheological evolution (t=2s), where the cloud map and contour lines represent the numerical values ​​of sediment volume fraction.

[0039] Figure 11 This is a simulation result of the spatial distribution of sediment volume fraction in flash floods and debris flows without considering rheological evolution (t=2s). The cloud map and contour lines represent the numerical values ​​of sediment volume fraction.

[0040] Figure 12 The figure shows the spatial distribution of rheological parameters of a flash flood and debris flow during the simulation (t=2s). Figure (A) shows the spatial distribution of yield stress, and Figure (B) shows the spatial distribution of viscosity. Detailed Implementation

[0041] The following examples further illustrate the method for simulating flash floods and debris flows based on rheological evolution characteristics provided by this invention. It should be noted that the following examples are only for further illustration of this invention and should not be construed as limiting the scope of protection of this invention. Any non-essential improvements and adjustments made by those skilled in the art based on the above-described invention will still fall within the scope of protection of this invention.

[0042] Example

[0043] This embodiment provides a detailed description of the flash flood and debris flow process simulation method based on rheological evolution characteristics described in this invention, as follows:

[0044] 1. The method for simulating flash floods and debris flows based on rheological evolution characteristics described in this invention is based on the relationship between the drag coefficient and Reynolds number of sediment particles carried by flash floods and debris flows during their movement (drag coefficient-Reynolds number relationship) derived in this invention. Here, the derivation process of the drag coefficient-Reynolds number relationship is explained first.

[0045] Since the mudslide slurry has Bingham fluid characteristics (see reference [1]: Huang Jiangcheng, Ou Guoqiang, Pan Huali, et al. Comparative study on rheological properties of typical soil slurry in mudslide-prone areas [J]. Journal of Sichuan University (Engineering Science Edition), 2015, 47(06): 49–53.), for Bingham fluid, its apparent viscosity η is expressed by equation (1):

[0046]

[0047] In equation (1), η is the apparent viscosity, η0 is the viscosity, and τ0 is the yield stress. denoted as shear rate.

[0048] Suppose a spherical particle with diameter d (hereinafter referred to as the particle) settles in an infinite Bingham fluid and reaches its maximum settling velocity v. max A dimensionless settling velocity ε* can be defined based on the units of the parameters in Bingham fluid, and its expression is shown in equation (2):

[0049]

[0050] In equation (2), ε* is the dimensionless settling velocity, η0 is the viscosity, τ0 is the yield stress, d is the particle diameter, and v max denoted as σ0, represents the maximum settling velocity of the particles in Bingham fluid.

[0051] Using the form of equation (1) and the dimensional analysis results of equation (2), the apparent viscosity is made dimensionless, resulting in the dimensionless apparent viscosity η. * Its expression is shown in equation (3):

[0052] η * =1+ε * (3)

[0053] In equation (3), η * For dimensionless apparent viscosity, ε * The value is the dimensionless settling velocity.

[0054] Based on dimensionless apparent viscosity η * The particle equivalent Reynolds number Re* is defined using equation (4):

[0055]

[0056] In equation (4), Re* ρ is the equivalent Reynolds number of the particle. f Let d be the density of the fluid containing the particle, d be the particle diameter, and v be the density of the fluid containing the particle. max η is the maximum settling velocity of particles in Bingham fluid, η0 is the viscosity, and ε * Let η be the dimensionless settling velocity, Re' be the Reynolds number of the particle, and η be the particle number. * It is a dimensionless apparent viscosity.

[0057] According to the research of Yue Xiang'an et al. (see reference [2]: Yue Xiang'an, Hao Jiangping, Chen Jialang. Resistance coefficient and settling velocity of solid particles in Bingham fluid [J]. Petroleum Drilling and Production Technology, 1993(1):1–8.), the resistance acting on particles in Bingham fluid is the equivalent Reynolds number Re of the particles. * The maximum settling velocity v of particles in Bingham fluid max The density ρ of the fluid containing the particles f And the particle diameter d and its related quantities are functions that can be expressed by equation (5):

[0058]

[0059] In equation (5), Re * ρ is the equivalent Reynolds number of the particle. f v is the density of the fluid containing the particle. max d is the maximum settling velocity of the particle in Bingham fluid, and d is the particle diameter.

[0060] Given that the particles have reached their maximum settling velocity in Bingham fluid, the drag they experience is the same as their own weight (see reference [3]: Zhang Xiaofeng, Liu Xingnian. River Dynamics [M]. Beijing: China Water Resources and Hydropower Press, 2010.), the expression for the drag coefficient of the particles is shown in equation (6):

[0061]

[0062] In equation (6), Cd is the drag coefficient of the particle, ρ p ρ is the particle density. f v is the density of the fluid containing the particle. max d is the maximum settling velocity of the particle in Bingham fluid, d is the particle diameter, and g is the acceleration due to gravity.

[0063] The experimental data of Valentik and Whitmore's sphere sedimentation experiment in Bingham fluid (see reference [4]: ​​Valentik L, Whitmore R L. The terminal velocity of spheres in Bingham plastics[J]. British Journal of Applied Physics, 1965, 16(8): 1197.) are as follows Figure 1 The data points are shown in the image. Figure 1 As shown, at lower Re * Below, Cd and Re * It exhibits a clear linear relationship on a log-log coordinate system, when Re * When Cd is large enough, it tends to a constant value, when Re * When the temperature is greater than 500, Cd is approximately 0.4. This trend is similar to the relationship observed in spherical particles settling in Newtonian fluids.

[0064] During the settling of particles in a Newtonian fluid, when the particle Reynolds number Re' is low, the relationship between the particle drag coefficient Cd and the particle Reynolds number Re' satisfies the relationship shown in equation (7):

[0065]

[0066] In equation (7), a and b are coefficients, Re' is the Reynolds number of the particle, and Cd is the drag coefficient of the particle.

[0067] from Figure 1 The experimental data shows that when Re * When the value is less than 10, the drag coefficient Cd of the particle is related to the equivalent Reynolds number Re of the particle. * It exhibits a significant linear relationship on a double logarithmic coordinate system. Therefore, by adopting the form of equation (7) and calibrating with the experimental results from the aforementioned reference [4], the stage (Re) is obtained. * <10) The relationship between the drag coefficient and the equivalent Reynolds number of the particles is shown in equation (8):

[0068]

[0069] In equation (8), Re' is the Reynolds number of the particle, and Cd is the drag coefficient of the particle.

[0070] Existing research indicates that when the settling velocity of spherical particles in Bingham fluid is high, the flow field around the particles becomes significantly asymmetric, leading to a negative wake behind them. In this case, the drag acting on the spherical particles mainly consists of pressure drag and turbulent friction drag. However, at lower Reynolds numbers, no significant separation occurs between the fluid and the spherical particles. A low Reynolds number indicates that viscous forces dominate the flow, resulting in a smooth laminar flow where the fluid adheres tightly to the surface of the spherical particles, thus reducing the likelihood of a separation zone forming.

[0071] The central angle of the low-pressure region projected behind the spherical particle is defined as the separation angle θ. It can be observed that for any spherical particle moving in the fluid, the separation angle ranges from 0 ≤ θ ≤ 2π. This is due to the particle's Reynolds number Re' (or the aforementioned equivalent Reynolds number Re). * The separation angle θ reflects the relative influence of turbulent drag and laminar viscous force. As the Reynolds number Re' of the particles increases, the influence of turbulent drag also increases. Therefore, the separation angle θ is directly related to the Reynolds number Re' of the particles. Existing studies have shown that the rate of increase of the separation angle θ decreases as the Reynolds number Re' of the particles increases (see reference [3]: Zhang Xiaofeng, Liu Xingnian. River Dynamics [M]. Beijing: China Water Resources and Hydropower Press, 2010.). Therefore, it is assumed that the separation angle θ has a certain influence on the equivalent Reynolds number Re' of the particles. * The derivative is inversely proportional to the equivalent Reynolds number of the particle, as shown in equation (9):

[0072]

[0073] In equation (9), C is a constant, and Re * The equivalent Reynolds number of the particle. This is the derivative of the separation angle with respect to the equivalent Reynolds number of the particle.

[0074] Re at the moment when the fluid and spherical particles are about to separate * Let Re *1 . When Re * Re *2 At this point, the vortex around the spherical particle fully develops, and a negative wake is formed. Therefore, the equation must satisfy the boundary conditions shown in equation (10):

[0075]

[0076] according to Figure 1 Experimental data analysis in the study determined Re *1 The value is 10, Re *2 The value is 500. Therefore, the boundary condition shown in equation (10) can be expressed by equation (11):

[0077]

[0078] In equation (11), θ is the separation angle, Re * is the equivalent Reynolds number of the particle.

[0079] The magnitude of turbulent resistance depends on the projected area of ​​the spherical particles in the separation zone in the flow direction (see reference [3]: Zhang Xiaofeng, Liu Xingnian. River Dynamics [M]. Beijing: China Water Resources and Hydropower Press, 2010.), and the area A affected by this turbulent resistance is... t As given by equation (12):

[0080]

[0081] In equation (12), A t Let d be the area affected by turbulent drag, and d be the particle diameter. Let be the projected area coefficient of the separation zone, and its expression is shown in equation (13):

[0082]

[0083] In equation (13), This is the projected area coefficient of the separation zone.

[0084] Once flow separation occurs, laminar viscous forces are limited to the surface region outside the separation zone. Therefore, the effective area A of the laminar viscous forces is... l As given by equation (14):

[0085]

[0086] In equation (14), A l θ is the effective area of ​​the laminar viscous force, θ is the separation angle, and d is the particle diameter.

[0087] Since the surface of the spherical particle is divided into regions where laminar viscous forces and turbulent drag are acting, the total drag acting on the particle is the sum of these two forces, as shown in equation (15):

[0088] F d =F t +F l (15)

[0089] In equation (15), F d F is the total resistance acting on the particle. t F is the turbulent resistance acting on the particle. l This refers to the laminar viscous force acting on the particles.

[0090] Therefore, the total resistance F acting on the particle d It can be expressed by equation (16):

[0091]

[0092] In equation (16), F d A is the total resistance acting on the particle. l A is the effective area of ​​laminar viscous force. t ρ is the area affected by turbulent drag. f v is the density of the fluid containing the particle. max Re represents the maximum settling velocity of particles in Bingham fluid. * is the equivalent Reynolds number of the particle.

[0093] Therefore, the semi-empirical formula for the drag coefficient of particles can be expressed by equation (17):

[0094]

[0095] In equation (17), Cd is the drag coefficient of the particle, and Re * Let θ be the equivalent Reynolds number of the particle, and θ be the separation angle. This is the projected area coefficient of the separation zone.

[0096] In the derivation of equation (17), no explicit introduction was made. Figure 1 Experimental data within the range of 10 ≤ Re* ≤ 500. However, the theoretical calculation results of equation (17) are plotted on... Figure 1 As can be seen from this, Figure 1 The curve in the figure is shown as the "Theoretical Derivation of Formula". Figure 1 The calculated values ​​(curves) are in high agreement with the experimental data (data points), thus verifying the applicability of Equation (17) in calculating the drag coefficient of particles in Bingham fluid.

[0097] 2. The rationality of the drag coefficient-equivalent Reynolds number relationship shown in equation (17) was verified by experiments.

[0098] To verify the rationality of the drag coefficient-equivalent Reynolds number relationship obtained by step 1 as shown in equation (17), experimental data on non-Newtonian fluid (mud) dam failure reported by Komatina and Jovanovic (see reference [5]: Komatina D, Jovanovic M. Experimental study of steady and unsteady free surfaceflows with water-clay mixtures[J]. Journal of Hydraulic Research, 1997, 35(5): 579–590.) were used for verification. The experimental setup and methods used for verification are described in reference [5], and the specific schematic diagram of the setup is shown in the figure. Figure 2 As shown.

[0099] To reduce computational load, the coarse-grained model proposed by Zhao et al. (see reference [6]: Zhao T, Utili S, Crosta G B. Rockslide and impulse wave modelling in the vajont reservoir by dem-cfd analyses[J]. Rock Mechanics and Rock Engineering, 2016, 49(6): 2437–2456.) was used to simulate an experiment with a volumetric sand content Cv = 20.1%, and the process diagram of the mud front position changing over time was plotted. The results are as follows: Figure 3 As shown. By Figure 3 It can be seen that the numerical simulation and experimental data are in good agreement, indicating that it is feasible to use the drag coefficient-equivalent Reynolds number relationship shown in Equation (17) obtained in step 1 above to calculate the motion problem of particles in non-Newtonian fluids by incorporating the CFD-DEM coupling method.

[0100] 3. Flowchart of the method for simulating flash floods and debris flows based on rheological evolution characteristics as described in this invention.

[0101] The flash flood and debris flow simulation method based on rheological evolution characteristics described in this invention is based on establishing a generalized numerical model of the target channel, and then using a CFD-DEM coupled model (CFDEM coupling engine) to simulate the transport process of sediment particles in the generalized numerical model under the action of water flow or debris flow slurry. The calculation flow of the method described in this invention is as follows: Figure 4 As shown, Figure 4 This demonstrates the computation process within one time step, as follows:

[0102] S1. The hilly gully to be simulated for the flash flood and debris flow process is taken as the target gully. The target gully includes interconnected debris flow formation areas and debris flow flow areas. The debris flow flow area gully has loose deposits. Based on the topographic information of the target gully and the gradation information of the sediment particles in the loose deposits, a generalized numerical model of the target gully is established in the open source software Gmsh.

[0103] S2, the CFD-DEM coupled model is used to simulate the transport process of sediment particles in a generalized numerical model under the action of water flow or debris flow slurry. The steps are as follows:

[0104] S21, Read the generalized numerical model of the target channel.

[0105] S22, the DEM module executes the DEM solution loop. In this loop, the particle neighbor list is first read, and contact and collision between particles are identified. Based on Hertzian contact theory and Mindlin and Deresiewicz theory, the nonlinear elastic normal and tangential forces between particles and between particles and boundaries are calculated, respectively. Then, the motion state of each particle is determined according to Newton's second law, and the position and velocity of each particle are updated. The position and velocity information of the particles are then passed to the coupling module of the CFD-DEM coupled model. The DEM solution loop process in this step is existing technology, and the open-source LIGGGHTS is used as the calculation software to simulate the movement of sediment particles, referring to existing technology.

[0106] S23: The particle position and velocity information obtained from the DEM module is transferred to the CFD module via the coupling module of the CFD-DEM coupled model. A CFD solution loop is executed in the CFD module. During this loop, the existing PISO algorithm is used to solve the fluid continuity equation and the Navier-Stokes equation, updating the flow field information. Once the flow field converges, the flow field information is transferred back to the coupling module of the CFD-DEM coupled model. The CFD solution loop process in this step is existing technology; referring to existing technology, the open-source OpenFOAM is used as the computational software to simulate water flow.

[0107] S24, in the coupling module of the CFD-DEM coupled model, the fluid velocity and mesh porosity at the location of each particle are obtained, the relative velocity of the particle fluid is calculated, and the density value obtained by weighted average of the fluid volume fraction of the mesh where the particle is located is used to determine whether the particle is located in the water flow. If the particle is located in the water flow, the rheological parameters are calculated according to the porosity of each mesh. The rheological parameters include the yield stress τ0 and the viscosity η0. If the particle is located in the air, the drag coefficient of the particle is calculated according to the Di Felice formula. For particles located in the water flow, it is determined whether the particle has started. If the particle has not started, the drag coefficient of the particle is calculated according to the Di Felice formula. If the particle has started, the equivalent Reynolds number Re of the particle is calculated according to equation (4). * Then, the separation angle and projected area coefficient are calculated using equations (11) and (13), and the drag coefficient of the particle is calculated using equation (17). Finally, the total drag F experienced by each particle is calculated. d The total resistance information experienced by the particles is transmitted back to the DEM module.

[0108]

[0109] In the formula, Re * ρ is the equivalent Reynolds number of the particle. f Let d be the density of the fluid containing the particle, d be the particle diameter, and v be the density of the fluid containing the particle. maxη is the maximum settling velocity of particles in Bingham fluid, η0 is the viscosity, and ε * Let η be the dimensionless settling velocity, Re' be the Reynolds number of the particle, and η be the particle number. * The apparent viscosity is dimensionless, and θ is the separation angle. is the projected area coefficient of the separation zone, and Cd is the drag coefficient of the particles.

[0110] This step is the main improvement and innovation of the present invention compared to the prior art.

[0111] S25, repeat the operations of S22 to S24, and continue the calculation for the next time step.

[0112] 4. The method described in this invention will be explained in detail below through specific examples.

[0113] (1) Establish a generalized numerical model of the target channel

[0114] Using a hilly gully for simulating a flash flood and debris flow process as the target gully, the target gully includes interconnected debris flow formation zones and debris flow circulation zones. The debris flow circulation zone contains loose sediment. Based on the topographic information of the target gully and the particle size distribution of the loose sediment, a generalized numerical model of the target gully is established in the open-source software Gmsh. In constructing the generalized numerical model, the location information of the midpoints, lines, surfaces, and volumes of the target gully are first compiled into an input file containing all necessary geometric information. This input file is then imported into the open-source software Gmsh, and a mesh is created using the mesh module according to the required model dimensions. The output is then saved to obtain the generalized numerical model of the target gully. Specifically, based on the topographic information of an actual hilly gully (target gully) and the particle size distribution of the loose sediment, a flash flood and debris flow moving bed model is established, including an upstream water tank, an upstream debris flow formation zone, and an interconnected debris flow circulation zone. A schematic diagram is shown below. Figure 5 As shown. The model uses an upstream water tank for water supply. When the fluid volume in the water tank is large enough, it can achieve an approximate steady-state flow at the outlet (see reference [7]: Nie Yipin, Lan Ling, Wang Xiekang. Numerical test on sediment deposition characteristics at the mouth of mountain torrents in high mountain canyons [J]. Engineering Science and Technology, 2024, 56(1): 237–244., reference [8]: Wang Youbiao. Study on impact force of debris flow on bridge piers [D]. Chengdu: Southwest Jiaotong University, 2019.). Based on the objective of this invention, namely, the simulation of mountain torrent debris flow process considering rheological evolution characteristics, the upstream water tank stores clear water. During the process of water flow eroding the bed surface, the rheological evolution is caused by the change in sediment content, and the overall flow characteristics gradually approach those of mountain torrent debris flow. That is, this step simulates the transport process of sediment particles in the generalized numerical model under the action of water flow using the CFD-DEM coupled model.

[0115] In the aforementioned flash flood and debris flow bed model, the gradient of both the debris flow formation zone and the debris flow circulation zone is 5%. The lengths of the upstream debris flow formation zone and the downstream debris flow circulation zone are 1.0 m and 3.0 m, respectively, and the width of both is 0.2 m. A 0.2 m thick layer of silt is laid in the debris flow circulation zone as an erodible bed surface, with a particle size of 20 mm (using a coarse-particle model, a particle size scaling ratio of 10, equivalent to 2 mm silt particles).

[0116] (2) Simulation of flash flood and debris flow process

[0117] To compare the impact of varying volumetric sediment content of debris flows on the rheological evolution of debris flows during simulation, an experimental and control group were set up for comparison, as detailed below:

[0118] Experimental group (simulation considering rheological evolution): After establishing the generalized numerical model of the target channel in step (1), the initial rheological parameters of the fluid in the upstream tank are set as follows: η0 = 10 -3 Pa·s, τ0=0Pa, the process of mountain torrent and debris flow is simulated according to the method of step S2 in step 3. That is, in the simulation of mountain torrent and debris flow, the rheological parameters η0 and τ0 are adjusted in real time according to the change of volumetric sand content Cv in different fluid grids. The force on the particles is calculated by the semi-empirical formula of the resistance coefficient as shown in formula (17) obtained in this invention, and the data of the position, velocity and other changes of sediment particles in the generalized numerical model with time are obtained.

[0119] Control group (simulation without considering rheological evolution): After establishing the generalized numerical model of the target channel in step (1), the rheological parameters are preset before simulating the flash flood and debris flow process (the initial parameters of the fluid in the upstream tank are set as: η0 = 10). -3 Pa·s, τ0=0Pa), regardless of the change in volumetric sediment concentration Cv during the simulation, the rheological parameters η0 and τ0 remain constant. Referring to step S2 in step 3 above, the flash flood and debris flow process is simulated. During the simulation, the Di Felice formula built into the open-source software is directly used to calculate the drag coefficient and total drag of the particles, without needing to first calculate the equivalent Reynolds number Re of the particles according to equation (4) as in the experimental group. * Then, the separation angle and projected area coefficient are calculated using equations (11) and (13), and the drag coefficient of the particles is calculated using equation (17). Finally, the data on the position, velocity, etc. of the sediment particles in the generalized numerical model are obtained by simulation over time.

[0120] (3) Comparison of simulation results

[0121] ① Characteristics of gully erosion during flash floods and debris flows

[0122] Figure 6 and Figure 7 This describes the evolution of gully erosion characteristics during the movement of flash floods and debris flows. Figure 6 These are simulation results considering the rheological evolution after a flash flood or debris flow comes into contact with particles. Figure 7 These are simulation results of flash floods and debris flows after they come into contact with sediment, without considering rheological evolution. Figure 6 It is known that, considering rheological evolution, debris flows, after entering an erodible bed, tend to move towards the bottom of the bed, forming deep scour pits at the leading edge of the bed. Due to the change in the volumetric sediment content (Cv) of the debris flow, which leads to rheological evolution, the energy dissipation caused by viscosity increases during the downstream movement of the debris flow and sediment mixture, thus slowing down the flow's velocity. And from... Figure 7 It can be seen that, in the simulation without considering rheological evolution, the interaction between the water and sediment phases remains unchanged, and the velocity of the debris flow is not affected by the additional viscosity caused by the change in the sediment concentration (Cv) of the debris flow. As the debris flow moves, thanks to the decrease in the relative velocity between the debris flow and the sediment, although the large volumetric sediment concentration (Cv) of the debris flow in some fluid domains causes rheological evolution, the influence of sediment particles on the movement of the debris flow is small. At t=5s, regardless of whether rheological evolution is considered, the debris flow has already contacted the bottom of the moving bed region.

[0123] ② Characteristics of changes in the elevation of the gully bed and water surface during flash floods and debris flows

[0124] Figure 8 and Figure 9 This describes the characteristics of gully bed and water surface elevation changes during flash floods and debris flows. Figure 8 It considers the rheological characteristics of the movement of mountain torrents and debris flows, including the changes in gully bed and water surface elevation. Figure 9 This refers to the characteristics of gully bed and water surface elevation changes in flash floods and debris flows without considering rheology. Figures 8-9 It can be seen that under the continuous scouring of flash floods and debris flows, scour pits of varying depths and widths were formed at the front end of the moving bed region. With the extension of the simulation time, regardless of whether rheological evolution was considered, the scour pits showed a trend of vertical deepening and widening along the flow path. For example... Figure 8 As shown, after the rheological evolution of sediment at the contact bed surface in flash floods and debris flows, the resistance to movement in the flow direction increases, making it more inclined to move towards the bottom of the moving bed, thus intensifying the erosion of lower sediment particles and significantly increasing the depth of the leading-edge scour pit. Figure 9As shown, in simulations that do not consider rheological evolution, the debris flow experiences relatively little resistance in the flow direction, forming scour pits that are long but shallow. Furthermore, the sediment stirred up by the scour forms a "sand mound"-like topographic feature on the bed surface. As bed scour intensifies, the high concentration of sediment carried by the debris flow in the "sand mound" area enhances the fluid viscosity. The reduced flow velocity of the debris flow promotes sediment deposition, which in turn reduces the sediment content of the debris flow, decreasing its viscosity and increasing its sediment-carrying capacity. This leads to a gradual increase in the height of the "sand mound" and its slow downstream migration, as shown in the example. Figure 8 As shown. Without considering rheological evolution, the movement of sediment is constrained by its inherent properties and its relative velocity with the debris flow. Therefore, when a debris flow comes into contact with the bed sediment, gravity causes a deflection of the flow direction, forming a scour pit and a subsequent "sand ridge" at the leading edge of the moving bed. Without considering rheological evolution, the debris flow crossing the "sand ridge" can still maintain a high velocity, carrying sediment downstream, causing the height of the "sand ridge" to undergo a process of first increasing and then decreasing, as shown in the diagram. Figure 9 As shown.

[0125] From the trend of water surface elevation characteristics (see...) Figures 8-9 (Water surface elevation data in the data), when the sediment content of flash floods and debris flows changes, whether or not rheological evolution is considered will significantly affect the water surface line simulation results. For example... Figure 8 As shown, in the simulation considering rheological evolution, the leading edge of the flash flood / debris flow entering the moving bed region carries a large amount of sediment, concentrating its energy and forming a deep scour pit near the inlet. At this point, most of the sediment particles carried by the flash flood / debris flow move along the bed surface. Subsequent flash floods / debris flows entering the moving bed region cover the water flow that has already entered the region, with lower sediment content and viscosity, enabling them to cross the "sand ridge" with less energy loss, resulting in a significant increase in water surface elevation. Figure 9 As shown, in simulations that do not consider rheological evolution, flash floods and debris flows create scour pits that are relatively long but shallow. Subsequent flash floods and debris flows encounter less obstruction along their path, making it difficult for them to achieve significant elevation changes, resulting in a relatively low water level. However, due to the limited obstruction along the path, the downstream propagation speed from the highest point of the water surface is faster than in simulations that consider rheological evolution.

[0126] The spatial distribution simulation results of sediment particle volume fraction in flash floods and debris flows at t=2s are as follows: Figure 10 and Figure 11 As shown, the cloud map and contour lines represent the numerical values ​​of sediment volume fraction. Figure 10 The simulation results take into account rheological evolution. Figure 11 These are simulation results that do not consider rheological evolution. For example... Figure 10As shown, in simulations considering rheological evolution, the sediment volume fraction is higher at the leading edge of flash floods and debris flows, while it is lower in the subsequent fluid immediately following the leading edge, resulting in sediment concentration at the leading edge. Figure 11 As shown, in the simulation without considering rheological evolution, since the rheological parameters do not change, less energy is consumed during the fluid flow due to internal viscosity. This results in the fluid front moving at a significantly higher velocity than the sediment particles, leading to a lower volume fraction of sediment particles in the fluid front and a higher volume fraction of sediment particles in the subsequent fluid.

[0127] Figure 12 The figures show the spatial distribution of rheological parameters of debris flow during the simulation of the debris flow process (t=2s). Figure (A) shows the spatial distribution of yield stress τ0, and Figure (B) shows the spatial distribution of viscosity η0. Since the drag coefficient-equivalent Reynolds number relationship shown in the derived equation (17) is closely related to the volumetric sediment content (Cv) of the debris flow, the two figures are similar in shape. In general, the rheological parameters show an increasing trend from top to bottom in the vertical direction. This is because under non-extreme hydrodynamic conditions, most sediment particles move along the bed surface, resulting in a higher sediment content in the debris flow near the bed surface, and a significant increase in the rheological parameter values. Figure 12 This study also revealed that the leading edge of flash floods and debris flows with a distance exceeding 2.5 meters along the course of the flow exhibits high yield stress and viscosity, indicating a high concentration of sediment particles in this area and the potential to cause significant damage to the surrounding environment and structures under high-speed movement. The volumetric sediment content (Cv) decreases in the main body of the flash flood and debris flow, and the yield stress and viscosity also decrease accordingly. However, this area has a relatively high water level, and the potential risk of inundation should be a key focus.

[0128] Based on the simulation results of flash floods and debris flows with or without considering rheological evolution, it can be seen that whether or not rheological evolution is considered in the simulation process will have significantly different effects on the simulation results, including the characteristics of gully erosion during the movement of flash floods and debris flows, and the characteristics of changes in the elevation of the gully and water surface during the movement of flash floods and debris flows.

[0129] Because the volumetric sediment content of a debris flow changes during its transport, leading to changes in its rheological parameters, the method described in this invention fully considers the impact of changes in volumetric sediment content on the rheological characteristics of debris flows during transport. By adjusting the rheological parameters of the debris flow in real time during the simulation, the changes in rheological parameters caused by changes in volumetric sediment content can be reflected in the simulation, thereby more realistically simulating the spatial distribution of sediment particles during the transport of debris flows.

Claims

1. A method for simulating flash floods and debris flows based on rheological evolution characteristics, characterized in that, Includes the following steps: S1. The hilly gully to be simulated for the flash flood and debris flow process is taken as the target gully. The target gully includes the interconnected debris flow formation area and debris flow flow area. The debris flow flow area gully has loose deposits. A generalized numerical model of the target gully is established based on the topographic information of the target gully and the gradation information of the sediment particles in the loose deposits. S2, the CFD-DEM coupled model is used to simulate the transport process of sediment particles in a generalized numerical model under the action of water flow or debris flow slurry. The steps are as follows: S21, Read the target channel generalized numerical model; S22, the DEM module executes the DEM solving loop. In the DEM solving loop, the particle neighbor list is first read and the contact and collision between particles are identified. The force on each particle is calculated. Then the motion state of each particle is determined, the position and velocity of each particle are updated, and the position and velocity information of the particles are passed into the coupling module of the CFD-DEM coupled model. S23, the position and velocity information of the particles obtained by the DEM module is passed to the CFD module through the coupling module of the CFD-DEM coupled model. The CFD solution loop is executed in the CFD module. In the CFD solution loop, the fluid control equation is solved and the flow field information is updated. When the flow field converges, the flow field information is passed to the coupling module of the CFD-DEM coupled model. S24. In the coupling module of the CFD-DEM coupled model, the fluid velocity and mesh porosity at the location of each particle are obtained, the relative velocity of the particle fluid is calculated, and it is determined whether the particle is in the water flow. If the particle is in the water flow, the rheological parameters including yield stress and viscosity are calculated according to the porosity of each mesh. If the particle is in the air, the drag coefficient of the particle is calculated according to the DiFelice formula. For particles in the water flow, it is determined whether the particle has started. If the particle has not started, the drag coefficient of the particle is calculated according to the DiFelice formula. If the particle has started, the equivalent Reynolds number of the particle is calculated according to equation (4). Then, the separation angle and projected area coefficient are calculated through equations (11) and (13), and the drag coefficient of the particle is calculated through equation (17). Finally, the total resistance experienced by each particle is calculated, and the total resistance information of the particles is transmitted back to the DEM module; In the formula, Re * ρ is the equivalent Reynolds number of the particle. f Let d be the density of the fluid containing the particle, d be the particle diameter, and v be the density of the fluid containing the particle. max η is the maximum settling velocity of particles in Bingham fluid, η0 is the viscosity, and ε * Let η be the dimensionless settling velocity, Re' be the Reynolds number of the particle, and η be the particle number. * The apparent viscosity is dimensionless, and θ is the separation angle. is the projected area coefficient of the separation zone, and Cd is the drag coefficient of the particle; S25, repeat the operations of S22 to S24, and continue the calculation for the next time step.

2. The method for simulating flash floods and debris flows based on rheological evolution characteristics according to claim 1, characterized in that, In step S2, when using the CFD-DEM coupled model to simulate the transport process of sediment particles in the generalized numerical model under the action of water flow, the water flow rate or debris flow slurry flow rate is determined by analyzing the historical rainfall data at the target channel.

3. The method for simulating flash floods and debris flows based on rheological evolution characteristics according to claim 1, characterized in that, In step S22, calculating the force on each particle refers to calculating the nonlinear elastic normal force and tangential force between particles and between particles and boundaries based on Hertzian contact theory and Mindlin and Deresiewicz theory, respectively; after determining the force on each particle, the motion state of each particle is determined according to Newton's second law.

4. The method for simulating flash floods and debris flows based on rheological evolution characteristics according to claim 1, characterized in that, In step S23, the fluid control equations are the fluid continuity equation and the NS equation. The PISO algorithm is used to solve the fluid control equations.

5. The method for simulating flash floods and debris flows based on rheological evolution characteristics according to any one of claims 1 to 4, characterized in that, In step S1, a generalized numerical model of the target gully is established in the open-source software Gmsh based on the topographic information of the target gully and the gradation information of the sediment particles in the loose deposits.

Images