Numerical simulation method of earth-rock dam disaster caused by cracks in the through dam culvert

Abstract

The application provides a numerical simulation method for a disaster-causing process of a dam-penetrating culvert pipe crack of an earth-rock dam, comprising: obtaining initial definite solution conditions of the dam-penetrating culvert pipe and the dam body, and a mapping relationship between a particle discrete classification of a soil body and physical parameters; performing non-steady seepage finite element calculation by using a current permeability coefficient at a current time step to obtain a centroid seepage slope of each unit; comparing and judging the centroid seepage slope with a critical seepage slope of a corresponding region according to a preset dam body region division; dynamically updating the unit permeability coefficient based on a judgment result, recording an adjustment number, and feeding back to a next time step for closed loop iteration; and determining a unit damage state and a leakage channel through condition based on the adjustment number. The application can automatically present a physical law of damage backtracking evolution without complex test calibration, and provides an evaluation basis for dam safety early warning.

Description

Technical Field

[0001] This invention relates to the field of water conservancy engineering safety evaluation technology, and more specifically, to a numerical simulation method for the disaster-causing process of cracks in culverts penetrating earth-rock dams. Background Technology

[0002] Culverts, as common water conveyance structures in earth-rock dams, are vulnerable to damage due to their penetrating nature. In practical engineering, culverts often develop cracks due to uneven foundation settlement or material aging. When high-pressure water seeps into the dam through these cracks, it triggers contact erosion at the pipe-soil interface and piping within the soil, leading to soil loss, changes in permeability, and the formation of seepage channels. Therefore, studying the mechanism and evolution of dam damage caused by culvert cracks and seepage is of mechanical research value for accurately assessing dam seepage stability and preventing engineering hazards.

[0003] Currently, numerical simulations of seepage disaster processes in earth-rock dams primarily employ continuous erosion rate models to characterize the loss of soil particles. These methods rely on complex pitting corrosion tests to calibrate empirical parameters such as erosion rates. Furthermore, when simulating the seepage field, they typically treat the seepage parameters across the entire computational domain as static constants or use a single critical hydraulic gradient threshold for global state determination. However, in actual seepage failure processes, the stress boundary conditions differ significantly between the soil within the dam body, the soil-pipe interface, and the downstream outlet. Moreover, the loss of local soil particles leads to nonlinear deterioration of the surrounding hydraulic conditions. Existing unified parameter models struggle to accurately capture the anisotropic, localized progressive failure characteristics at these weak points, resulting in limited accuracy in predicting the seepage deformation propagation rate and spatial evolution path under complex conditions.

[0004] In summary, existing methods struggle to simultaneously address the coupled response mechanisms of local medium property dynamics evolution and macroscopic seepage field redistribution when dealing with progressive seepage failure under complex boundary conditions. Therefore, there is an urgent need to investigate a numerical simulation method capable of characterizing the self-accelerating effect of local damage and anisotropic failure features to improve the reliability of predicting the evolution of dam seepage disasters under complex scenarios. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide an improved technical solution for the numerical simulation method of the disaster-causing process of cracks in the culverts penetrating earth-rock dams.

[0006] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: A numerical simulation method for the disaster-causing process of cracks in culverts penetrating earth-rock dams includes: Obtain the initial boundary conditions of the culvert and dam body, as well as the pre-established mapping relationship between the discrete classification of soil particles and physical property parameters; Based on the initial boundary conditions and mapping relationship, unsteady seepage finite element calculation is performed using the current permeability coefficient of each calculation unit at the current time step to obtain the seepage head of each calculation node and the centroid seepage gradient of each calculation unit. Based on the pre-defined dam area division, the centroid seepage gradient of each calculation unit is compared with the critical seepage gradient of the corresponding area to obtain the discrimination result; Based on the discrimination results and mapping relationship, the current permeability coefficient of the calculation unit that meets the triggering condition is dynamically updated, and the number of times the permeability coefficient of the calculation unit is adjusted is recorded. The updated current permeability coefficient is fed back to the unsteady seepage finite element calculation of the next time step for closed-loop iteration. Based on the number of permeability coefficient adjustments, the damage status of each calculation unit inside the dam and the connectivity of seepage channels are determined, and the simulation results of the disaster-causing process are output.

[0007] In an exemplary embodiment, for the pipe-soil contact area, the centroidal seepage gradient of each calculation unit is compared with the critical seepage gradient of the corresponding area to obtain a discrimination result, specifically including: The centroidal seepage gradient of the calculation unit located in the pipe-soil contact area is vector-decomposed into tangential and normal components along the pipe wall. The tangential component was compared with the critical gradient of the contact surface scour, and the normal component was compared with the critical gradient of the internal piping. When the tangential component is greater than the critical gradient of the contact surface scour, or the normal component is greater than the critical gradient of the internal piping, the judgment result is determined to meet the triggering condition.

[0008] In one exemplary embodiment, establishing the mapping relationship between particle size and critical permeability gradient further includes: When the dam area is divided into the downstream outlet zone, the slope angle between the downstream dam slope and the horizontal plane is introduced as a force balance parameter to calculate the critical slope drop that takes into account the tendency of soil particles to slide down the dam slope.

[0009] In one exemplary embodiment, during the closed-loop iteration process, the step size of the time step is dynamically adjusted based on the changes in the dam body's phreatic surface obtained from the unsteady seepage finite element calculation. In the initial stage of the closed-loop iteration, the first step length is used; when the seepage head variation between adjacent time steps is less than the preset head stability threshold, it is determined that the dam body seepage surface tends to be stable, and the time step length is increased to a second step length value that is greater than the first step length value.

[0010] In one exemplary embodiment, the simulation results of the disaster-causing process are output, including: Based on the seepage head, centroid seepage gradient, current permeability coefficient, and damage state of each calculation node obtained from the unsteady seepage finite element calculation at each time step, an isohead distribution map, a seepage gradient contour map, and a dam damage distribution evolution map are generated and output to provide dam safety early warning information.

[0011] The beneficial effects of this invention are: this invention does not require complex experimental calibration, and can automatically present the physical laws of the retrospective evolution of damage, providing an assessment basis for dam safety early warning. Attached Figure Description

[0012] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with this application and, together with the description, serve to explain the principles of this application.

[0013] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, for those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0014] Figure 1 This is a flowchart illustrating the steps of a numerical simulation method for the disaster-causing process of cracks in the culvert penetrating an earth-rock dam, as described in an embodiment of this application.

[0015] Figure 2 This is a flowchart illustrating the steps involved in offline construction of a physical evolution parameter library to drive subsequent dynamic simulations, as described in an embodiment of this application.

[0016] Figure 3 This is a flowchart illustrating the steps involved in performing finite element calculations according to an embodiment of this application.

[0017] Figure 4 This is a flowchart illustrating the steps of the positive feedback iteration mechanism in an embodiment of this application.

[0018] Figure 5 This is a flowchart illustrating the steps for determining the damage status of each computational unit inside the dam body in an embodiment of this application. Detailed Implementation

[0019] To enable those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present application, and not all embodiments. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative effort should fall within the scope of protection of the present application.

[0020] Furthermore, in the embodiments of this application, the words "exemplary," "for example," etc., are used to indicate that they are examples, illustrations, or descriptions. Any embodiment or design scheme described as "exemplary" in the embodiments of this application should not be construed as being better or more advantageous than other embodiments or design schemes. Specifically, the use of the term "exemplary" is intended to present the concept in a concrete manner.

[0021] References to "one embodiment" or "some embodiments" as described in this specification mean that one or more embodiments of this application include a specific feature, structure, or characteristic described in connection with that embodiment. Therefore, the phrases "in one embodiment," "in some embodiments," "in other embodiments," "in still other embodiments," etc., appearing in different parts of this specification do not necessarily refer to the same embodiment, but rather mean "one or more, but not all, embodiments," unless otherwise specifically emphasized. The terms "comprising," "including," "having," and variations thereof mean "including but not limited to," unless otherwise specifically emphasized.

[0022] In this application embodiment, "at least one" refers to one or more, and "more than one" refers to two or more. "And / or" describes the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A alone, A and B simultaneously, and B alone, where A and B can be singular or plural. The character " / " generally indicates that the preceding and following related objects are in an "or" relationship. "At least one of the following" or similar expressions refer to any combination of these items, including any combination of single or plural items. For example, at least one of a, b, or c can represent: a, b, c, ab, ac, bc, or abc, where a, b, and c can be single or multiple.

[0023] Example 1: This embodiment provides a numerical simulation method for the disaster-causing process of cracks in culverts penetrating earth-rock dams, combined with... Figure 1 As shown, the method includes the following steps 101-105: Step 101: Obtain the initial boundary conditions of the culvert and dam body, as well as the pre-established mapping relationship between the discrete classification of soil particles and physical property parameters; determine the initial permeability coefficient of each calculation unit based on the mapping relationship, and use it as the current permeability coefficient.

[0024] Specifically, initial boundary conditions are used to provide the basic boundaries and initial states for subsequent unsteady seepage calculations. In actual engineering modeling, initial boundary conditions are implemented by setting various physical boundaries. For example, the cracks in the culverts penetrating the dam are set as equal-head boundaries, the upstream reservoir water level is set as a constant-head boundary, the downstream dam slope is set as a free seepage section, and the permeable layer boundary at the bottom of the dam is set as a flow surface with no flow interaction.

[0025] Furthermore, the pre-established mapping relationship between soil particle discrete gradation and physical property parameters is specifically manifested as a parameter evolution library generated offline. This mapping relationship establishes a deterministic mathematical correlation between particle size distribution and soil porosity, critical seepage gradient, and permeability coefficient by discretizing fine soil particles according to particle size, providing a basis for determining the values ​​of permeability parameters during subsequent dynamic simulations.

[0026] Step 102: Based on the initial boundary conditions and mapping relationship, perform unsteady seepage finite element calculation using the current permeability coefficient of each calculation unit at the current time step to obtain the seepage head of each calculation node and the centroid seepage gradient of each calculation unit.

[0027] In this embodiment, the unsteady seepage finite element calculation is a discrete solution process that progresses over time. Initially, each calculation unit uses the permeability coefficient of the soil before particle loss. Upon entering any current time step, the program constructs a global permeability matrix using the current permeability coefficients of each calculation unit, and solves the seepage differential equation using the initial boundary conditions to obtain the head distribution field of the entire dam's computational domain. Based on this, the scalar or vector permeability gradient at the centroid of each calculation unit is extracted through spatial derivative operations. The centroidal permeability gradient characterizes the local hydraulic gradient of that calculation unit under hydraulic action at the current time step and is a mechanical parameter driving seepage deformation in the soil.

[0028] Step 103: Based on the preset dam body area division, compare and judge the centroid seepage gradient of each calculation unit with the critical seepage gradient of the corresponding area to obtain the judgment result.

[0029] Accordingly, due to the spatial differences in stress state and failure mechanism of the soil surrounding the culvert, this method adopts a zonal discrimination logic. The preset dam body region division specifically includes the internal soil region, the pipe-soil contact surface region, and the downstream outlet region. After obtaining the centroid seepage gradient, the program determines the spatial region to which the current calculation unit belongs and calls the seepage resistance criterion specific to that region. By performing numerical comparison calculations, it determines whether the centroid seepage gradient reaches or exceeds the critical seepage gradient of the corresponding region, which serves as the physical trigger condition for whether soil particles inside the unit will begin to flow out, and then outputs a Boolean type discrimination result.

[0030] Step 104: Based on the discrimination results and mapping relationship, dynamically update the current permeability coefficient of the calculation unit that meets the triggering condition, record the number of times the permeability coefficient of the calculation unit is adjusted, and feed the updated current permeability coefficient back to the unsteady seepage finite element calculation of the next time step for closed-loop iteration.

[0031] In some optional implementations, the trigger condition is met when the determination result indicates that the centroidal permeability gradient of a certain calculation unit is greater than the critical permeability gradient. At this time, the program, based on the mapping relationship, jumps the permeability coefficient of the calculation unit from the current level to a higher value corresponding to the next level of particle loss, and simultaneously increments the permeability coefficient adjustment counter configured for that unit by 1. If the centroidal permeability gradient is not greater than the critical permeability gradient, the permeability coefficient remains unchanged.

[0032] Furthermore, the updated permeability coefficient is reassembled into the overall permeability matrix of the next time step, constructing a two-way positive feedback closed loop of seepage field calculation and soil damage evolution. That is, local particle loss leads to an increase in permeability coefficient, and the increase in permeability coefficient causes the redistribution of seepage field, which may accelerate the deterioration of hydraulic conditions of surrounding units, reflecting the physical law of spontaneous expansion of seepage failure.

[0033] Step 105: Based on the number of permeability coefficient adjustments, determine the damage status of each calculation unit inside the dam body and the connectivity of the seepage channels. When it is determined that the seepage channels are connected, terminate the closed-loop iteration and output the simulation results of the disaster-causing process.

[0034] Specifically, the number of times the permeability coefficient is adjusted assesses the degree of loss of fine particles inside the computing unit and is a quantitative indicator characterizing the accumulation of local damage. By comparing the number of adjustments with a preset threshold, the computing unit can be accurately classified into different damage stages, such as safe, near-destruction, or complete destruction.

[0035] Furthermore, by tracking the spatial topological connectivity of computational units in a completely destroyed state, it can be determined whether the seepage channels are connected along the pipe-soil contact surface at the macroscopic level. Finally, the computational data from each time step are used to generate simulation results such as isohead distribution maps, seepage gradient contour maps, and dam damage distribution evolution maps, providing engineering data support for dam safety early warning and hazard assessment.

[0036] Example 2: This embodiment provides a method for constructing a soil particle discrete grading and physical property parameter mapping model. Using conventional geotechnical test data, an offline physical evolution parameter library is constructed to drive subsequent dynamic simulations. Figure 2 As shown, the construction method includes the following steps 201-2022: Step 201: The pre-established mapping relationship between discrete soil particle classification and physical property parameters is obtained by obtaining the initial physical parameters of the soil, which include at least the median particle size and the initial porosity.

[0037] Specifically, the median particle size can be obtained by conducting particle size distribution tests on the dam construction soil and plotting the particle size distribution curve. For the initial porosity, conventional soil mechanics methods are used to calculate the initial porosity value based on the dry density of the soil sample, the specific gravity of the soil particles, and the density of water. This calculation process is a standard technique in this field and will not be elaborated upon here.

[0038] Step 202: Divide the fine particles in the soil that are smaller than the median particle size into multiple particle grades in order of increasing particle size; determine the representative particle size and cumulative loss mass fraction corresponding to each particle grade as the data basis for constructing the mapping relationship.

[0039] In this embodiment, the discrete hierarchical mapping relationship is constructed based on the following physical assumptions: it is assumed that only particle loss occurs per unit volume of soil under seepage without considering the replenishment of new particles; it is also assumed that the particle loss process occurs sequentially in order of increasing particle size, that is, the loss of larger particle sizes will only be triggered after the loss of smaller particle sizes is completed. Based on the above assumptions, fine particles below the median particle size are divided into m particle levels, and the representative particle size of the i-th level and its corresponding cumulative loss mass fraction are determined.

[0040] Step 2021 establishes the mapping relationship between particle size and permeability coefficient, specifically including: for each particle size, based on the initial porosity and the cumulative loss mass fraction corresponding to that particle size, and according to the principle of mass conservation, calculate the current porosity after all the loss of that particle size.

[0041] Specifically, the initial solid volume fraction is obtained by calculating the difference between the value 1 and the initial porosity; the remaining mass ratio is obtained by calculating the difference between the value 1 and the cumulative loss mass fraction corresponding to that particle size; the current solid volume fraction is obtained by calculating the product of the initial solid volume fraction and the remaining mass ratio; and the current porosity after all the particles of that particle size have been lost is obtained by calculating the difference between the value 1 and the current solid volume fraction, as shown in the following formula: ;

[0042] In the formula: Indicates the first The current porosity of the soil after all particles have been lost; 1 represents the overall unit volume or unit proportion; Indicates the initial porosity of the soil; Determined by the particle size distribution curve, its value is equal to the percentage of the mass of soil particles with a particle size less than or equal to the representative particle size of the i-th grade, relative to the mass of all fine particles below the median particle size.

[0043] Step 2022: Based on the current porosity, the representative particle size corresponding to the particle size, and the physical property parameters of water, calculate the permeability coefficient corresponding to the particle size after loss according to the capillary model and Darcy's law, so as to establish the mapping relationship between each particle size and the permeability coefficient.

[0044] Specifically, the numerator is calculated by multiplying the square of the current porosity, the square of the particle size, and the specific gravity of water; the denominator is calculated by multiplying the constant 18, the dynamic viscosity coefficient of water, and the difference between the current solid volume fraction (i.e., the value of 1) and the current porosity; the ratio of the numerator to the denominator is then calculated to obtain the permeability coefficient for the corresponding particle size distribution. The formula is as follows: ;

[0045] In the formula: Indicates the first The permeability coefficient corresponding to the loss of particles; Indicates the first The current porosity of the soil after all particles have been lost; Indicates the first Representative particle size of the particles; Indicates the specific gravity of water; 18 represents the constant multiplier derived from the capillary model of parallel equal-diameter circular tubes; This represents the dynamic viscosity coefficient of water; 1 represents the total unit volume.

[0046] As an alternative, considering that the permeability coefficients calculated based on the above formula do not reflect the complexity of the actual soil pore structure, their calculated values ​​may be too high. In scenarios where experimental conditions permit, physical test data can be used to calibrate and correct the theoretically calculated permeability coefficients at various levels, thereby optimizing the accuracy of the parameter mapping library.

[0047] Optionally, a mapping relationship between particle size and critical seepage gradient is established, specifically including: for each particle size, based on the current porosity and the representative particle size corresponding to that particle size, the critical seepage gradient corresponding to that particle size in each dam area division is calculated respectively. Specifically, based on the structural characteristics and stress states of different parts of the dam, the dam body is divided into the internal soil zone, the soil-pipe interface zone, and the downstream outlet zone. Since different zones correspond to different failure mechanisms, this embodiment employs differentiated critical seepage gradient calculation logics for each of the three zones.

[0048] When the dam area is divided into an internal soil zone, the unit volume cohesion of the soil in the initial physical parameters is introduced as a force balance parameter to calculate the corresponding critical gradient for internal piping.

[0049] Specifically, calculate the current solid volume fraction, which is the product of the current porosity, the effective unit weight of soil particles, and the tangent of the internal friction angle of the soil, minus the value of 1. Add this product to the cohesion per unit volume to obtain the resistance numerator. Calculate the ratio of the resistance numerator to the unit weight of water to obtain the critical gradient for internal piping. The formula is as follows: ;

[0050] In the formula: Indicates the first After the loss of particles, the critical seepage gradient of internal piping in the internal soil zone. 1 represents the total unit volume; Indicates the current porosity; Indicates the effective unit weight of soil particles; Indicates the internal friction angle of the soil; This represents the cohesion per unit volume of soil, and its value is equal to the cohesion measured by the direct shear test. (Unit: kPa) divided by 1m (i.e.) This is to ensure that the dimensions are consistent with other terms in the formula (kN / m³); It indicates the density of water.

[0051] Optionally, when the dam body area is divided into a soil-pipe contact area, the crack width, which increases dynamically with the increase of the representative particle size, is introduced as a stress balance parameter to calculate the corresponding critical slope for scour of the contact surface.

[0052] Specifically, the product of constant 4, the current solid volume fraction, the representative particle size, the effective unit weight of soil particles, and the tangent of the internal friction angle of the soil is used as the numerator of the resistance term; the product of constant 3, the crack width, and the unit weight of water is used as the denominator of the dynamic term; the ratio of the resistance numerator to the dynamic denominator is calculated to obtain the critical slope of the contact surface scour, and the calculation formula is as follows: ;

[0053] In the formula: Indicates the first After the loss of particles, the critical seepage gradient of the contact surface corresponding to the pipe-soil contact area is determined. 4 represents the constant multiplier generated from the force equilibrium derivation; 1 represents the total unit volume; Indicates the current porosity; This indicates the particle size; Indicates the effective unit weight of soil particles; Indicates the internal friction angle of the soil; 3 represents the constant multiplier generated during the derivation of force equilibrium; It represents the width of the micro-cracks at the contact surface between the culvert and the soil. It dynamically increases with the progression of particle loss level, and qualitatively characterizes the evolution trend of the dynamic expansion of micro-channels at the pipe-soil contact surface due to the loss of fine particles.

[0054] In this embodiment, the mapping relationship between the crack width and the representative grain size is set as follows: ;

[0055] In the formula: Indicates the first The width of the cracks corresponding to the loss of particles; This represents the crack propagation ratio coefficient, which is set to 3 in this embodiment. Those skilled in the art can make adaptive adjustments based on engineering conditions such as culvert material, contact surface roughness, and relative displacement between the pipe and the soil. Indicates the first The representative particle size of each particle size class. The above linear mapping relationship reflects the physical understanding that when particles of a certain size are completely lost, the microscopic channels available for water flow at the contact surface approximately expand to the same order of magnitude as the particle size of that size. It indicates the density of water.

[0056] Optionally, establishing the mapping relationship between particle size and critical seepage gradient also includes: when the dam area is divided into the downstream outlet zone, introducing the slope angle between the downstream dam slope and the horizontal plane as a force balance parameter to calculate the critical slope gradient considering the tendency of soil particles to slide down the dam slope.

[0057] Specifically, the product of the cosine of the slope angle and the tangent of the internal friction angle is calculated, and the sine of the slope angle is subtracted to obtain the slope resistance coefficient; the sum of the products of the cosine of the slope angle, the sine of the slope angle, and the tangent of the internal friction angle is calculated to obtain the slope dynamic coefficient; the product of constant 4, the current solid volume fraction, the representative particle size, the effective unit weight of soil particles, and the slope resistance coefficient is used as the numerator; the product of constant 3, the fracture width, the unit weight of water, and the slope dynamic coefficient is used as the denominator; the ratio of the numerator to the denominator is calculated to obtain the critical slope gradient for slope slide, and the calculation formula is as follows: ;

[0058] In the formula: Indicates the first After the loss of primary particles, the critical infiltration gradient of the slope corresponding to the downstream outlet area. 4 represents the constant multiplier generated from the force equilibrium derivation; 1 represents the total unit volume; Indicates the current porosity; This indicates the particle size; Indicates the effective unit weight of soil particles; This indicates the angle between the downstream dam slope and the horizontal plane; Indicates the internal friction angle of the soil; 3 represents the constant multiplier generated from the force equilibrium derivation; This represents the width of the microcrack, and its value is determined according to the mapping relationship above. Sure; It indicates the density of water.

[0059] Example 3: This embodiment provides a strategy for preparatory work and dynamic step size adjustment in unsteady seepage calculations. This strategy refines the execution details of finite element calculations, combining... Figure 3 As shown, the steps include 301-302 as follows: Step 301, before performing the unsteady seepage finite element calculation, also includes: dividing the dam body into finite element meshes to obtain multiple calculation units; among which, the pipe-soil contact area and the downstream outlet area where the culvert is in contact with the dam body are subjected to mesh densification processing.

[0060] Specifically, this embodiment involves the setting of the underlying basic model for finite element calculation. The entire seepage field solution domain is spatially discretized using conventional triangular elements to generate multiple computational elements, and the head function within each element is approximated using linear interpolation. Subsequent unsteady seepage calculations utilize conventional iterative methods to solve the overall finite element equations. The aforementioned mesh generation and equation solving process belongs to existing finite element processing methods in this field.

[0061] Based on this, considering that the pipe-soil contact area and the downstream outlet area are high-incidence areas of seepage deformation and structurally weak surfaces, and that their local hydraulic gradient changes have high-frequency characteristics, mesh refinement was applied to these two areas. By increasing the density of computational cells per unit area, higher-resolution local head field and seepage gradient data were obtained, providing data support for subsequent damage state assessment.

[0062] Step 302: During the closed-loop iteration process, the step size of the time step is dynamically adjusted according to the changes in the dam body saturation surface obtained by the unsteady seepage finite element calculation. In the initial stage of the closed-loop iteration, the first step size is used. When the seepage head variation between adjacent time steps is less than the preset head stability threshold, it is determined that the dam body saturation surface tends to be stable, and the time step is increased to a second step size greater than the first step size.

[0063] In this embodiment, the time step represents the characteristic sampling interval for the redistribution of the seepage field. Within each time step, the permeability coefficient of each computational unit undergoes at most one transition, representing the completion of a certain level of particle loss event when the hydraulic conditions continuously meet a certain critical gradient during that period. The selection of the time step must be less than the characteristic time for wetting surface adjustment and ensure a smooth transition of seepage field data between adjacent time steps.

[0064] To achieve dynamic adjustment of the step size, a unified process for determining the state of the wetting surface is established: The seepage head amplitude data and seepage flow amplitude data for two adjacent time periods are acquired. This seepage head amplitude data is compared with a preset head stability threshold. The head stability threshold and flow stability threshold can be determined by those skilled in the art based on the actual engineering accuracy requirements. For example, the head stability threshold can be set to the order of 0.1% of the maximum head value, and the flow stability threshold can be set to the order of 1% of the initial seepage flow. The seepage flow amplitude data is then compared with a preset flow stability threshold. Based on the comparison results, the following parallel branch logic is executed: Branch 1: When it is determined that the dam body's seepage surface has not stabilized (e.g., the seepage head amplitude data is greater than or equal to the head stabilization threshold, or the seepage flow amplitude data is greater than or equal to the flow stability threshold), it is determined that it is in the initial stage of closed-loop iteration, and a smaller first step length value is adopted (e.g., set to 0.5~1 time units). This setting is used to capture the changes in the seepage field and seepage surface in the initial stage at a higher sampling frequency.

[0065] Branch 2: When it is determined that the dam's seepage surface is stabilizing (e.g., the seepage head amplitude is less than the head stabilization threshold, and the seepage flow amplitude is less than the flow rate stabilization threshold), the time step is increased to a second step value greater than the first step value (e.g., set to 2-4 time units, with the upper limit not exceeding the set tolerance limit). This setting is used to reduce redundant iterative calculations and improve the computational efficiency of the entire numerical simulation process.

[0066] Example 4: As an optional implementation, this embodiment simulates the anisotropic weak surface of the dam body by introducing a three-zone differential discrimination mechanism and vector decomposition logic of the pipe-soil contact surface, and constructs a positive feedback iterative mechanism for local damage-driven seepage field redistribution.

[0067] Combination Figure 4 As shown, the mechanism includes: a pre-defined dam body area division comprising an internal soil area, a pipe-soil contact surface area, and a downstream outlet area; for the internal soil area, the critical gradient of internal piping determined based on the internal piping mechanism is used as the critical permeability gradient for the corresponding area; for the pipe-soil contact surface area, the critical gradient of contact surface scour determined based on the contact surface scour mechanism is used as the tangential critical permeability gradient for the corresponding area, and the critical gradient of internal piping determined based on the internal piping mechanism is used as the normal critical permeability gradient for the corresponding area; for the downstream outlet area, the critical gradient of slope slippage determined based on the slope slippage failure mechanism is used as the critical permeability gradient for the corresponding area.

[0068] In practical implementation, the dam body zoning is defined based on physical spatial boundaries. For example, the area from the outer wall of the culvert penetrating the dam to the predetermined channel expansion instability thickness (the reference value for this thickness is usually in the range of 0.5m to 1.0m) is designated as the pipe-soil contact surface area; the local boundary area at a predetermined horizontal distance from the downstream dam face (such as 0.1m to 0.5m) is designated as the downstream outlet area; and the remaining area inside the dam body, excluding the pipe-soil contact surface area and the downstream outlet area, is designated as the internal soil area. Through the above spatial zoning, spatial positioning benchmarks are provided for the occurrence of different failure mechanisms, and the corresponding critical seepage gradient values ​​are invoked according to the mapping relationship.

[0069] Optionally, when comparing and judging the centroidal permeability gradient of each calculation unit with the critical permeability gradient of the corresponding region, the following steps are included: determining the damage reduction coefficient of adjacent units based on the number of times the permeability coefficient of adjacent units is adjusted; dynamically correcting the critical permeability gradient of the corresponding region using the damage reduction coefficient of adjacent units to obtain the corrected critical permeability gradient; and comparing and judging the centroidal permeability gradient of each calculation unit with the corrected critical permeability gradient.

[0070] This embodiment constructs a local state positive feedback mechanism in the closed loop of seepage field evolution. Specifically, it extracts the number of times the permeability coefficient of adjacent cells within the topological neighborhood of the target computational cell is adjusted, and calculates the damage reduction coefficient of adjacent cells based on the number of adjustments.

[0071] Specifically, the number of permeability coefficient adjustments for each adjacent cell within the topological neighborhood of the target computational cell is extracted, and the maximum value is selected as the maximum number of adjustments in the neighborhood. Based on the ratio of the maximum number of adjustments in the neighborhood to the complete destruction threshold, the damage reduction coefficient of the adjacent cells is calculated. In this embodiment, the formula for calculating the damage reduction coefficient of the adjacent cells is as follows: ;

[0072] In the formula: This represents the damage reduction coefficient between adjacent units, and its value ranges from greater than 0 to no greater than 1. 1 represents the baseline value corresponding to the undamaged state; This represents the damage transfer strength coefficient, which is set to 0.8 in this embodiment. Those skilled in the art can make adaptive adjustments based on the erosion resistance characteristics of the soil in actual engineering projects. This represents the maximum number of times the permeability coefficient has been adjusted among all adjacent cells within the topological neighborhood of the target computational cell; This indicates the threshold for complete destruction.

[0073] When the calculated result of the above formula is less than the preset lower limit (0.2 in this embodiment), the damage reduction coefficient of adjacent units is set to this preset lower limit to avoid the critical seepage gradient being corrected to an excessively small non-physical value. The product of the critical seepage gradient of the corresponding region and the damage reduction coefficient of adjacent units is calculated to obtain the corrected critical seepage gradient. This process simulates the physical phenomenon of local particle loss leading to a weakening of the erosion resistance of the neighborhood.

[0074] When performing the above comparison and discrimination, a unified regional attribute determination action is adopted for the calculation units of different regions, and the corresponding parallel comparison branch logic is executed according to the determination result: Branch 1: When the calculation unit is determined to be located in the internal soil zone or the downstream outlet zone, extract the centroid seepage gradient scalar of the calculation unit and compare the scalar with the modified critical seepage gradient of the corresponding area.

[0075] Branch 2: When the calculation unit is located in the pipe-soil contact area, the weak surface in this area exhibits anisotropic characteristics. Perform the following steps: Optionally, for the pipe-soil contact area, the centroidal seepage gradient of each calculation unit is compared with the critical seepage gradient of the corresponding area to obtain the discrimination result. Specifically, this includes: vector decomposing the centroidal seepage gradient of the calculation unit located in the pipe-soil contact area into tangential and normal components along the pipe wall; comparing the tangential component with the critical gradient of the contact surface scour, and comparing the normal component with the critical gradient of internal piping; when the tangential component is greater than the critical gradient of the contact surface scour, or the normal component is greater than the critical gradient of internal piping, the discrimination result is determined to meet the triggering condition.

[0076] In the above vector decomposition process, the tangential component represents the scouring driving force of water flow along the outer wall of the culvert, and the normal component represents the piping driving force of water flow perpendicular to the outer wall of the culvert and penetrating into the dam body.

[0077] Specifically, for each calculation unit located within the pipe-soil contact area, the local normal direction of the culvert's outer wall at that unit is determined based on the geometric relationship between its centroid position and the culvert axis. The centroidal seepage gradient vector is projected onto the normal direction to obtain the normal component; the centroidal seepage gradient vector is projected onto a direction perpendicular to the normal direction and along the pipe wall to obtain the tangential component. In comparisons, the critical seepage gradient is used as the benchmark. This directional discrimination mechanism overcomes the limitations of single scalar discrimination and improves the sensitivity of trigger condition identification at the contact surface.

[0078] After obtaining the discrimination results of each computing unit within the entire computing domain, a uniform penetration coefficient update process is performed on all computing units: Optionally, the current permeability coefficient of the calculation unit that meets the triggering condition is dynamically updated, including: when the centroidal permeability gradient is greater than the critical permeability gradient of the corresponding region, the triggering condition is determined to be met; according to the mapping relationship, the current permeability coefficient of the calculation unit that meets the triggering condition is updated to the permeability coefficient corresponding to the next level of particle loss; when the centroidal permeability gradient is not greater than the critical permeability gradient of the corresponding region, the current permeability coefficient of the calculation unit is kept unchanged.

[0079] Specifically, for the calculation unit executing branch two, its centroidal permeability gradient is greater than the critical permeability gradient, meaning that either the tangential or normal component exceeds its limit. For calculation units that meet the triggering condition, its current permeability coefficient k is... i The transition is updated to the next-level permeability coefficient k in the mapping relationship. i+1 Simultaneously, the permeability coefficient adjustment count counter for that unit is incremented by 1. After completing the update actions for all computational units within the current time step, the updated permeability coefficient set is fed back to the computation program to assemble the overall permeability matrix for the next time step, driving the recalculation of the seepage head distribution and completing the closed-loop iterative process.

[0080] Example 5: This embodiment provides a quantitative characterization of dam damage and a disaster determination mechanism based on the number of permeability coefficient adjustments, and gives the multi-dimensional result output format required for engineering early warning.

[0081] Combination Figure 5 As shown, optionally, the damage state of each calculation unit inside the dam body is determined based on the number of permeability coefficient adjustments, including: determining the critical failure threshold and the complete failure threshold based on the total number of particle grades in the discrete grading of soil particles; comparing the number of permeability coefficient adjustments of the calculation unit with the critical failure threshold and the complete failure threshold respectively; if the number of permeability coefficient adjustments is greater than the critical failure threshold but less than the complete failure threshold, the damage state of the calculation unit is determined to be critical failure; if the number of permeability coefficient adjustments reaches the complete failure threshold, the damage state of the calculation unit is determined to be complete failure.

[0082] Specifically, the total number of particle classifications represents the number of grades into which fine particles below the median particle size are divided. The threshold calculation logic is as follows: calculate the sum of the total number of particle classifications and the value 1 to obtain the complete destruction threshold; calculate the product of the complete destruction threshold and a preset proportion coefficient (e.g., set to 60%) to obtain the near-destruction threshold.

[0083] For example, in a simplified normalized example, assuming the total number of particle gradations is set to 5, the calculated threshold for complete destruction is 6 times; the product of 6 and 60% yields 3.6 times, and in the actual judgment logic, the threshold for near-destruction is determined to be 4 times through a conventional rounding up operation. Through the above threshold comparison logic, the complex continuous physical damage process is discretized into a clear numerical state range, providing a discrete state basis for subsequent disaster judgment.

[0084] As an optional implementation, determining the damage state of each calculation unit inside the dam body and the continuity of the seepage channel also includes: when the damage state of the calculation unit is determined to be a completely destroyed state, setting the current permeability coefficient of the calculation unit to a preset maximum permeability coefficient value; monitoring the spatial distribution of calculation units in a completely destroyed state; when the calculation units in a completely destroyed state form a continuous chain of destroyed units along the pipe-soil contact surface of the culvert through the dam, from the crack in the culvert through the dam to the downstream outlet, determining that the seepage channel is connected, and terminating the closed-loop iteration.

[0085] Furthermore, when the cumulative simulation time of the closed-loop iteration reaches the preset maximum simulation duration, and no through-seepage channel has formed inside the dam body, the dam body is determined to be in a seepage-safe state during that simulation period, and the closed-loop iteration is terminated. The maximum simulation duration can be set according to the time scale of interest in the actual project, and those skilled in the art can make adaptive adjustments as needed.

[0086] During the discrimination and update process at each time step, the damage state threshold determination is performed first: when the number of times the permeability coefficient of the computing unit is adjusted reaches the complete destruction threshold, the stepwise permeability coefficient transition update based on the mapping relationship is no longer performed, but the processing logic of the complete destruction state is entered.

[0087] In this embodiment, when a computational unit reaches a state of complete failure, it indicates that all soil particles with a median diameter or smaller within the unit have been lost, resulting in substantial damage to the soil skeleton. At this point, the permeability coefficient corresponding to the unit after the loss of the last stage of particles is obtained, multiplied by a preset magnification factor (e.g., 10), to obtain a preset maximum permeability coefficient value, which is then assigned to the current permeability coefficient of the unit. This assignment operation aims to numerically characterize the high water conductivity of the large pore channels formed after the skeleton is destroyed.

[0088] Furthermore, a conventional spatial topology tracing algorithm is used to monitor the connectivity of computational units in a completely destroyed state. When it is determined that the connected domain composed of completely destroyed units forms a continuous topological path connecting the upstream crack and the downstream outlet in space, the seepage channel is considered to be connected. Channel connection indicates a sudden change in the hydraulic conditions inside the dam body, and the seepage field no longer conforms to the linear seepage premise of Darcy's law. Therefore, the system triggers a termination command, ending the current unsteady seepage finite element closed-loop iteration process.

[0089] Optionally, the simulation results of the disaster-causing process can be output, including: the seepage head, centroid seepage gradient, current permeability coefficient and damage state of each calculation node obtained from the unsteady seepage finite element calculation at each time step, and the generation and output of isohead distribution map, seepage gradient contour map and dam damage distribution evolution map, so as to provide dam safety early warning information.

[0090] Specifically, at each time step of the closed-loop iteration and after the simulation terminates, the numerical solution results of each computational unit are extracted. Using spatial interpolation and contour line generation algorithms, the seepage head data of each computational node are processed to generate contour maps of equal head lines; the centroid seepage gradient data of each computational unit are processed to generate contour maps of seepage gradient; and a damage distribution evolution map of the dam body is generated based on the number of times the permeability coefficient of each computational unit is adjusted.

[0091] Furthermore, to provide multi-dimensional safety early warning information, this embodiment also includes an optimized supplementary solution: selecting multiple spatial fixed points at preset intervals along the axis of the culvert penetrating the dam, extracting seepage head data at each time step for these fixed points, and generating a fixed-point seepage head process line; simultaneously, performing time integration or transient extraction on the seepage flow passing through the downstream outlet boundary to generate a seepage flow process line. The above multi-dimensional output intuitively reveals the spatial distribution and temporal evolution of dam damage, providing an assessment basis for engineers to predict the location of potential hazards, evaluate the rate of damage propagation, and formulate emergency response plans.

[0092] Those skilled in the art will understand that the specific output results of the above numerical simulation method (such as the seepage channel penetration time and damage propagation rate) are closely related to specific engineering conditions such as dam material parameters, culvert crack size, and reservoir water level. In practical applications, engineers can input the corresponding parameters for specific engineering scenarios and run the above closed-loop iterative process to obtain the corresponding disaster-causing process simulation results.

[0093] This scheme constructs a theoretical mapping relationship between discrete soil particle classification and physical property parameters, which can generate a physical evolution parameter library with only conventional geotechnical test data, thus lowering the threshold for engineering application of the model and improving the problem that existing numerical simulations rely on complex pitting corrosion test calibration parameters.

[0094] Meanwhile, the introduction of a three-zone differentiated critical slope criterion and the vector decomposition of the centroid slope at the pipe-soil contact surface help identify the failure priority order of different stress boundaries and characterize the anisotropic failure characteristics of weak surfaces, thus improving the difficulty of traditional models in depicting anisotropic failure and regional mechanism differences. Addressing the difficulty in coupling local damage with the macroscopic seepage field, a closed-loop iterative mechanism for dynamically updating the damage reduction coefficient and permeability coefficient of adjacent units is introduced. This mechanism can reflect the positive feedback effect of local erosion deterioration driving seepage field redistribution and can reproduce the physical evolution process of failure regression.

[0095] Furthermore, to address the lack of evaluation indicators, this invention uses the number of permeability coefficient adjustments as a benchmark to determine the discretized damage threshold and tracks the topological connectivity of completely destroyed units, thereby achieving quantitative characterization of the spatial distribution of damage inside the dam body and automated determination of the seepage channel connectivity status, providing data support for dam hazard early warning.

[0096] The sequence numbers of the above embodiments of the present invention are for descriptive purposes only and do not represent the superiority or inferiority of the embodiments.

[0097] In the above embodiments of the present invention, the descriptions of each embodiment have different focuses. For parts not described in detail in a certain embodiment, please refer to the relevant descriptions of other embodiments.

[0098] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details of the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A numerical simulation method for the disaster-causing process of cracks in culverts penetrating earth-rock dams, characterized in that, include: Obtain the initial boundary conditions of the culvert and dam body, as well as the pre-established mapping relationship between the discrete classification of soil particles and physical property parameters; Based on the initial boundary conditions and mapping relationship, unsteady seepage finite element calculation is performed using the current permeability coefficient of each calculation unit at the current time step to obtain the seepage head of each calculation node and the centroid seepage gradient of each calculation unit. Based on the pre-defined dam body area division, the centroid seepage gradient of each calculation unit is compared with the critical seepage gradient of the corresponding area to obtain the discrimination result; Based on the discrimination results and mapping relationship, the current permeability coefficient of the calculation unit that meets the triggering condition is dynamically updated, and the number of times the permeability coefficient of the calculation unit is adjusted is recorded. The updated current permeability coefficient is fed back to the unsteady seepage finite element calculation of the next time step for closed-loop iteration. Based on the number of permeability coefficient adjustments, the damage status of each calculation unit inside the dam body and the connectivity of seepage channels are determined, and the simulation results of the disaster-causing process are output.

2. The method according to claim 1, characterized in that, The pre-defined dam body area division includes the internal soil area, the pipe-soil contact area, and the downstream outlet area; For the internal soil area, the critical gradient of internal piping determined based on the internal piping mechanism is used as the critical permeability gradient of the corresponding area. For the soil-pipe contact area, the critical slope of the contact surface scour determined based on the contact surface scour mechanism is used as the tangential critical permeability slope of the corresponding area, and the critical slope of the internal piping determined based on the internal piping mechanism is used as the normal critical permeability slope of the corresponding area. For the downstream outlet area, the critical slope gradient determined based on the slope slip failure mechanism is used as the critical infiltration slope gradient for the corresponding area.

3. The method according to claim 1, characterized in that, When comparing the centroidal seepage gradient of each calculation unit with the critical seepage gradient of the corresponding region, the following steps are taken: The damage reduction factor of adjacent units is determined based on the number of adjustments made according to the permeability coefficient of adjacent units. The critical permeability gradient of the corresponding region is dynamically corrected by using the damage reduction coefficient of adjacent units, and the corrected critical permeability gradient is obtained. The centroidal seepage gradient of each calculation unit is compared with the modified critical seepage gradient for discrimination.

4. The method according to claim 1, characterized in that, Dynamically update the current permeability coefficient of the computational unit that meets the triggering conditions, including: When the centroidal permeability gradient is greater than the critical permeability gradient of the corresponding region, the triggering condition is determined to be met. Based on the mapping relationship, the current permeability coefficient of the calculation unit that meets the triggering condition is updated to the permeability coefficient corresponding to the loss of particles in the next stage; When the centroidal permeability gradient is not greater than the critical permeability gradient of the corresponding region, the current permeability coefficient of the calculation unit remains unchanged.

5. The method according to claim 1, characterized in that, Based on the number of permeability coefficient adjustments, the damage state of each computational unit inside the dam body is determined, including: Based on the total number of particle classifications in the discrete classification of soil particles, the thresholds for impending failure and complete failure are determined. The number of times the permeability coefficient of the calculation unit is adjusted is compared with the thresholds for near destruction and complete destruction, respectively; If the number of permeability coefficient adjustments is greater than the critical damage threshold but less than the complete damage threshold, then the damage state of the computing unit is determined to be critical damage. If the number of times the permeability coefficient is adjusted reaches the complete destruction threshold, the damage state of the computing unit is determined to be a complete destruction state.

6. The method according to claim 5, characterized in that, Determining the damage status of each calculation unit within the dam body and the connectivity of seepage channels also includes: When the damage state of the computing unit is determined to be complete destruction, the current permeability coefficient of the computing unit is set to the preset maximum permeability coefficient value. Monitor the spatial distribution of computing units in a completely destroyed state; When a computational unit in a completely destroyed state forms a continuous chain of destroyed units along the soil-to-pipe contact surface of the culvert, extending from the crack in the culvert to the downstream outlet, the seepage channel is determined to be connected, and the closed-loop iteration is terminated.

7. The method according to claim 1, characterized in that, The pre-established mapping relationship between discrete soil particle size distribution and physical property parameters is obtained through the following method: Obtain the initial physical parameters of the soil, which should include at least the median particle size and initial porosity. Fine particles in the soil that are smaller than the median particle size are classified into multiple particle sizes in order of increasing particle size. Determine the representative particle size and cumulative loss mass fraction corresponding to each particle size class as the data basis for constructing the mapping relationship.

8. The method according to claim 7, characterized in that, Establishing a mapping relationship between particle size distribution and permeability coefficient, specifically including: For each particle size, based on the initial porosity and the cumulative loss mass fraction corresponding to that particle size, the current porosity after all loss of that particle size is calculated according to the principle of mass conservation. Based on the current porosity, the representative particle size corresponding to this particle size, and the physical properties of water, the permeability coefficient corresponding to the particle size after loss is calculated according to the capillary model and Darcy's law, so as to establish the mapping relationship between each particle size and the permeability coefficient.

9. The method according to claim 8, characterized in that, Establishing a mapping relationship between particle size distribution and critical permeability gradient, specifically including: For each particle size class, based on the current porosity and the representative particle size corresponding to that particle size class, the critical seepage gradient corresponding to that particle size class in each dam area division is calculated respectively. When the dam area is divided into an internal soil area, the unit volume cohesion of the soil is introduced as a force balance parameter to calculate the corresponding critical gradient for internal piping. When the dam body area is divided into the soil-pipe contact area, the crack width, which increases dynamically with the increase of the representative particle size, is introduced as a stress balance parameter to calculate the corresponding critical slope for scour of the contact surface.

10. The method according to claim 1, characterized in that, Before performing the finite element calculation for unsteady seepage, the following steps are also included: The dam body is meshed using the finite element method to obtain multiple computational elements; Among them, the pipe-soil contact area and the downstream outlet area where the culvert is in contact with the dam body are subjected to grid densification treatment.

Images