A numerical analysis method for fluid-structure interaction under free drainage conditions of unsaturated slopes

Abstract

This invention discloses a numerical analysis method for fluid-structure interaction (FSI) of unsaturated slopes under free drainage conditions, belonging to the field of numerical analysis technology for slope engineering. The method includes: constructing an unsaturated multi-source parameter model of the slope region; building a dynamic boundary model with hysteresis response and gradual release characteristics based on boundary element water content, water potential gradient, and rainfall input information; jointly updating the internal water content gradient, local structural changes, and permeability coefficient field of the slope; correcting the seepage direction through a path enhancement factor; and identifying local preferential flow channels and seepage concentration areas. Iterative solutions are then obtained through adaptive time step adjustment, residual convergence control, and anomaly element correction strategies. This method improves the physical consistency of the free drainage boundary description, enhances the ability to identify preferential flow channels and local seepage concentration areas, and improves the computational stability and engineering applicability of FSI analysis of unsaturated slopes.

Description

Technical Field

[0001] This invention relates to the field of numerical analysis technology for slope engineering, specifically to a fluid-structure interaction numerical analysis method for unsaturated slopes under free drainage conditions. Background Technology

[0002] The fluid-structure interaction numerical analysis method for unsaturated slopes under free drainage conditions refers to a numerical calculation method for unsaturated soil slopes. Considering the dynamic changes in free drainage boundary conditions, it achieves a unified analysis of the internal water migration and mechanical response processes of the slope by coupling and co-solving the seepage field and the mechanical field. This method quantitatively assesses the stability of slopes under complex hydrological conditions by constructing an unsaturated multi-source parameter model, establishing a dynamic boundary drainage mechanism related to water state and time, adaptively modeling the seepage path, and combining fluid-structure interaction response calculation and solution control strategies.

[0003] In existing technologies, most analysis methods for unsaturated slopes are based on simplified boundary conditions or static parameter assumptions. For example, in boundary treatment, the slope surface is often simplified to a constant flux boundary or a fixed head boundary, which fails to reflect the drainage lag and gradual release characteristics during free drainage after rainfall stops. In seepage modeling, the permeability coefficient is often only related to the water content, without considering the impact of soil structure changes (such as deformation and crack development) on the seepage path and permeability, resulting in an inability to accurately describe the formation and evolution of local preferential flow channels. In fluid-structure interaction calculations, some methods adopt weak coupling or step-by-step calculation strategies, which suffer from problems such as coupling update lag and insufficient numerical stability, making it difficult to meet the analysis needs under strongly nonlinear conditions.

[0004] In practical engineering, such as the gradual transition of a slope surface from saturated to unsaturated after rainfall, drainage capacity often exhibits a gradual change from weak to strong. Locally, seepage concentration channels may form due to structural loosening or crack development, leading to local stress redistribution and deformation concentration. Existing methods struggle to simultaneously account for the dynamic coupling effects of these multiple factors, resulting in discrepancies between analytical results and actual working conditions.

[0005] Therefore, it is necessary to provide a numerical analysis method that can comprehensively consider the dynamic evolution of free drainage boundaries, adaptive adjustment of seepage paths, and collaborative updating of fluid-structure interaction responses, so as to improve the accuracy and engineering applicability of unsaturated slope stability analysis. Summary of the Invention

[0006] To address the above issues and overcome the shortcomings of existing technologies, this invention provides a fluid-structure interaction numerical analysis method for unsaturated slopes under free drainage conditions. This method includes the following steps:

[0007] Step S1: Modeling unsaturated multi-source parameters;

[0008] Step S2: Evolution of the dynamic boundary of free drainage;

[0009] Step S3: Penetration path modeling;

[0010] Step S4: Fluid-structure interaction response modeling;

[0011] Step S5: Fluid-structure interaction solution control.

[0012] Further, in step S1, the unsaturated multi-source parameter modeling is used to construct the basic parameter system of the unsaturated slope. Specifically, it involves unifying and structuring the soil physical parameters, water content parameters, and initial stress state of the slope area, combining field test data, expressing the soil water content, matrix suction relationship, and permeability characteristics, and performing parameter consistency correction and spatial mapping to obtain unsaturated parameter model data for subsequent calculations.

[0013] The unsaturated parameter model data specifically includes: basic soil physical property parameters, water content and matrix suction characteristic curve data, unsaturated permeability coefficient function data, and initial stress and initial water content data;

[0014] The basic physical property parameters of the soil are used to characterize the soil structure and compaction.

[0015] The moisture content versus matrix suction characteristic curve data are used to characterize the unsaturated moisture retention properties;

[0016] The unsaturated permeability coefficient function data is used to characterize the seepage capacity under different water content conditions;

[0017] The initial stress and initial water-bearing state data are used as initial condition inputs for fluid-structure interaction analysis.

[0018] Further, in step S2, the free drainage dynamic boundary evolution is used to characterize the dynamic change process of the boundary drainage capacity of the slope under free drainage conditions. Specifically, based on the water content of the surface area of ​​the slope, rainfall input information and water potential gradient distribution, the drainage capacity of the boundary area is dynamically determined, and a drainage response adjustment function related to time and local water content is constructed to continuously update the boundary drainage behavior and obtain dynamic boundary flux data.

[0019] The evolution of the free drainage dynamic boundary specifically adopts a free drainage dynamic boundary modeling method based on hysteresis response adjustment. In some preferred embodiments, a drainage response adjustment function can also be constructed to nonlinearly modulate the boundary drainage flux to reflect the hysteresis characteristics and gradual release characteristics in the drainage process.

[0020] The free drainage dynamic boundary modeling method based on hysteresis response adjustment includes the following steps:

[0021] Step S21: Boundary cell identification, used to determine the boundary regions participating in free drainage and establish the initial state. Specifically, the surface grid cells of the slope are screened to identify the boundary cells that meet the free drainage conditions, and the water content data, water potential data and rainfall input data corresponding to each boundary cell are extracted to obtain the initial boundary state data.

[0022] Step S22: Basic drainage flux calculation, used to obtain the basic drainage capacity of the boundary area under unsaturated conditions. Specifically, based on the unsaturated permeability coefficient function and the boundary water potential gradient, the basic drainage flux of each boundary unit is calculated and processed to obtain basic drainage flux data.

[0023] Step S23: Construction of drainage response regulation function, used to characterize the hysteresis and gradual release characteristics of drainage capacity during free drainage. Specifically, based on the water content data and time evolution information of the boundary unit, a nonlinear drainage response regulation function related to water content and time is constructed to modulate the drainage response process and obtain drainage response regulation function data.

[0024] Step S24: Dynamic boundary flux correction, used to dynamically adjust the basic drainage flux. Specifically, based on the drainage response adjustment function data, the basic drainage flux is nonlinearly corrected, and the flux is smoothly updated in combination with time continuity constraints to obtain the drainage flux distribution data of the boundary unit.

[0025] Step S25: Drainage capacity evolution update, used to construct the time change process of drainage capacity. Specifically, based on the time step advancement mechanism, the drainage flux of each boundary unit is continuously updated and processed to form the sequence data of drainage flux changing with time, and the drainage capacity changes with time data is obtained.

[0026] Step S26: Boundary water state response modeling, which is used to reflect the feedback effect of the drainage process on the boundary water state. Specifically, based on the change of drainage flux of the boundary unit, the water content of the boundary region is updated and calculated, and the trend of water state change over time is recorded to obtain the water state response data of the boundary region.

[0027] Step S27: Update the hysteresis adjustment parameters to improve the adaptive capability of the drainage response model. Specifically, based on historical response data during the drainage process, the key parameters in the drainage response adjustment function are dynamically adjusted to obtain drainage hysteresis adjustment parameter data.

[0028] The dynamic boundary flux data specifically includes boundary unit drainage flux distribution data, drainage capacity variation data over time, boundary region water content response data, and drainage hysteresis adjustment parameter data.

[0029] Further, in step S3, the seepage path modeling is used to describe the evolution process of the seepage path inside the slope under the dynamic boundary drive. Specifically, based on the dynamic boundary flux data obtained in step S2, the water content distribution and its spatial gradient inside the slope are analyzed. Combined with the local deformation state or structural change information of the soil, the permeability coefficient field is adaptively updated, and the seepage path is directionally corrected and continuously reconstructed to obtain seepage path evolution data.

[0030] The permeability path modeling specifically introduces a path enhancement factor based on the water cut gradient to weight and adjust the local seepage direction, including the following steps:

[0031] Step S31: Internal state field extraction, used to obtain basic state information related to the evolution of seepage path inside the slope and determine potential seepage concentration areas. Specifically, based on the dynamic boundary flux data obtained in step S2, the water content distribution inside the slope is updated and calculated, and the spatial gradient analysis of water content is performed. At the same time, the local structural change index of each unit is extracted, and the areas that meet the water content gradient and structural change threshold conditions are screened to obtain seepage candidate area data.

[0032] Step S32: Adaptive update of the permeability coefficient field, used to jointly adjust the permeability capacity according to the change in water content and the local structural response. Specifically, based on the unsaturated permeability coefficient function, the permeability coefficient of each unit is updated, and the water content gradient response term and the structural change response term are introduced to correct the permeability coefficient, so as to obtain the updated data of the spatial distribution of permeability coefficient.

[0033] Step S33: Seepage direction correction, which is used to enhance the guiding effect of water content gradient on seepage path direction. Specifically, a path enhancement factor is constructed based on the water content gradient of each unit, and the path enhancement factor is weighted and fused with the original seepage direction to correct the seepage direction and obtain seepage path direction distribution data.

[0034] Step S34: Local priority flow channel identification, used to identify continuous channel structures where seepage is concentrated within the slope. Specifically, based on the updated spatial distribution data of permeability coefficient and the distribution data of seepage path direction, the continuity of seepage between adjacent units is determined, and units that meet the conditions of enhanced permeability and directional consistency are connected and merged to obtain local priority flow channel identification data.

[0035] Step S35: Seepage concentration area identification, used to determine seepage concentration areas and construct the continuous evolution results of seepage paths. Specifically, based on the local priority flow channel identification data, the degree of seepage concentration in the area surrounding the channel is analyzed, and areas that meet the concentration conditions are identified. At the same time, the seepage paths are continuously reconstructed by combining time series information to obtain seepage concentration area identification data.

[0036] Step S36: Output the permeability path evolution results, which is used to structure and organize the results of each sub-step and form a unified output. Specifically, it involves associating and integrating the spatial distribution update data of the permeability coefficient, the distribution data of the permeation path direction, the identification data of local preferential flow channels, and the identification data of the permeation concentration area to obtain the permeability path evolution data.

[0037] The seepage path evolution data specifically includes: permeability coefficient spatial distribution update data, seepage path direction distribution data, local preferential flow channel identification data, and seepage concentration area identification data;

[0038] The spatial distribution update data of the permeability coefficient is used to characterize the permeability capacity distribution of each unit inside the slope at the current moment;

[0039] The seepage path direction distribution data is used to characterize the dominant seepage direction of each unit after direction correction;

[0040] The local preferential flow channel identification data is used to characterize the high-permeability continuous channel structure formed inside the slope;

[0041] The data identifying concentrated seepage areas is used to characterize key areas where local seepage accumulation is significant and may induce subsequent response changes.

[0042] Furthermore, in step S4, the fluid-structure interaction response modeling is used to perform coupled modeling of the mechanical response of the slope soil based on the seepage field update. Specifically, based on the seepage path evolution data obtained in step S3, the distribution of pore water pressure inside the soil is updated, and the stress field is calculated and processed in combination with the strength characteristics of unsaturated soil to realize the linkage update between the seepage field and the stress field and obtain fluid-structure interaction response data.

[0043] The fluid-structure interaction response data specifically includes: pore water pressure distribution data, soil stress and strain response data, local deformation and displacement field data, and strength parameter update data.

[0044] Furthermore, in step S5, the fluid-structure interaction solution control is used to perform numerical solution and computational control on the fluid-structure interaction process. Specifically, based on the fluid-structure interaction response data obtained in step S4, the time step, iteration convergence conditions and error control strategy in the calculation process are adjusted and processed, and multi-step iterative solution is performed to obtain slope stability analysis result data.

[0045] The slope stability analysis results data specifically include: slope stability evaluation index data, safety factor calculation results data, potential slip surface distribution data, and instability risk area identification data.

[0046] The beneficial effects achieved by the present invention using the above solution are as follows:

[0047] (1) In view of the technical problems in the existing unsaturated slope fluid-structure interaction analysis method, there are strong discreteness of multi-source parameters and the large deviation between the calculation results and the actual working conditions caused by the simplified treatment of boundary conditions. This scheme creatively adopts a modeling mechanism that combines unsaturated multi-source parameter modeling with free drainage dynamic boundary evolution. By uniformly and structurally expressing the soil physical parameters, water-bearing characteristic parameters and initial state, and introducing a dynamic boundary flux response function based on water-bearing state and time evolution, the continuous updating and nonlinear modulation of the boundary drainage capacity are realized, thereby improving the physical consistency and calculation accuracy of the boundary condition description.

[0048] (2) In view of the technical problem that existing unsaturated seepage analysis methods have fixed permeability coefficient fields or rely only on water content as a single variable, making it difficult to reflect local structural changes and the formation process of preferential flow channels, this scheme creatively adopts a permeability coefficient adaptive update model based on water content gradient and structural change response, and combines path enhancement factor to correct seepage direction. At the same time, it identifies and continuously reconstructs local preferential flow channels, and can identify key risk areas such as concentrated seepage zone at the toe of the slope or fracture flow path.

[0049] (3) In view of the technical problems in the existing fluid-structure interaction numerical analysis methods, such as the lag in the coupling update of the seepage field and stress field, the insufficient stability of the solution process, and the weak engineering applicability, this scheme creatively adopts a fluid-structure interaction response modeling and adaptive solution control mechanism driven by the evolution of seepage path. By dynamically updating the pore water pressure, stress, strain and strength parameters, and introducing the time step adaptive adjustment, residual convergence control and abnormal element correction strategy, the stable and efficient solution of the fluid-structure interaction process is realized, and the reliability and applicability of engineering calculation are improved. Attached Figure Description

[0050] Figure 1 A flowchart illustrating a fluid-structure interaction numerical analysis method for unsaturated slopes under free drainage conditions provided by this invention;

[0051] Figure 2 This is a flowchart illustrating the evolution of the free drainage dynamic boundary in step S2.

[0052] Figure 3 A flowchart illustrating the penetration path modeling for step S3.

[0053] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used together with the embodiments of the invention to explain the invention and do not constitute a limitation thereof. Detailed Implementation

[0054] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0055] In the description of this invention, it should be understood that the terms "upper", "lower", "front", "rear", "left", "right", "top", "bottom", "inner", "outer", etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0056] Example 1, see Figure 1 This invention provides a fluid-structure interaction numerical analysis method for unsaturated slopes under free drainage conditions, the method comprising the following steps:

[0057] Step S1: Modeling unsaturated multi-source parameters;

[0058] Step S2: Evolution of the dynamic boundary of free drainage;

[0059] Step S3: Penetration path modeling;

[0060] Step S4: Fluid-structure interaction response modeling;

[0061] Step S5: Fluid-structure interaction solution control.

[0062] By performing the above operations, this solution addresses the technical problems in existing unsaturated slope fluid-structure interaction analysis methods, such as the strong discreteness of multi-source parameters and the large deviation between calculation results and actual working conditions caused by simplified boundary conditions. This solution creatively adopts a modeling mechanism that combines unsaturated multi-source parameter modeling with free drainage dynamic boundary evolution. By uniformly and structurally expressing soil physical parameters, water-bearing characteristic parameters, and initial state, and introducing a dynamic boundary flux response function based on water-bearing state and time evolution, the continuous updating and nonlinear modulation of boundary drainage capacity can be achieved.

[0063] For example, in the free drainage process after rainfall, traditional methods usually simplify the boundary to a fixed drainage or zero flux boundary, which cannot reflect the actual law of drainage lag in the early stage and enhancement in the later stage. However, this scheme uses a drainage response adjustment function to characterize the drainage lag and gradual release process, so that the boundary flux gradually transitions from a restricted state to a stable release state, thereby improving the physical consistency and calculation accuracy of the boundary condition description.

[0064] Example 2, see Figure 1 This embodiment is based on the above embodiment. In step S1, the unsaturated multi-source parameter modeling is used to construct the basic parameter system of the unsaturated slope. Specifically, it involves unifying and structuring the soil physical parameters, water content parameters, and initial stress state of the slope area, combining field test data, expressing the soil water content, matrix suction relationship, and permeability characteristics, and performing parameter consistency correction and spatial mapping to obtain unsaturated parameter model data for subsequent calculations.

[0065] First, the slope area was divided into zones. Based on geological survey data, borehole data, and field sampling results, different soil layers were numbered and labeled, and corresponding soil unit datasets were established. For each soil unit, its basic physical parameters were collected, including natural density, dry density, void ratio, liquid limit, plastic limit, and particle size distribution parameters. Natural density was determined by on-site ring sampling, dry density was obtained by the drying method, void ratio was calculated based on density relationships, liquid limit and plastic limit were obtained through laboratory liquid and plastic limit tests, and particle size distribution parameters were obtained through sieve analysis or laser particle size analysis. For missing data, statistical averages of similar soil types were used for completion.

[0066] Secondly, regarding the construction of water content characteristic parameters, unsaturated soil-water characteristic tests were conducted on each soil unit. The pressure plate apparatus or filter paper method was preferably used to obtain matrix suction data under different water content conditions, and a correlation curve between water content and matrix suction was constructed based on the test data. During data fitting, the Van Genuchten model was preferably used for parameter expression, where relevant parameters could be solved by inversion using the least squares method. When test data was lacking for some soil units, typical empirical parameter ranges were selected based on the soil type and assigned values, and spatial supplementation was achieved through interpolation of neighboring units.

[0067] Furthermore, regarding the construction of unsaturated permeability characteristic parameters, the unsaturated permeability coefficient can be expressed as a function based on the above-mentioned relationship between water content and matrix suction, combined with the test results of saturated permeability coefficient. Among them, the saturated permeability coefficient can be obtained through indoor constant head or variable head permeability tests. For the permeability in the unsaturated range, it can be derived and calculated based on the soil-water characteristic curve, and the results can be smoothed to avoid local abrupt changes from causing unstable effects on subsequent calculations.

[0068] Regarding the construction of the initial stress state and initial water content state, the initial stress field distribution can be calculated using a layered gravity loading method based on the slope geometry and soil self-weight conditions. The initial water content state can be determined through on-site moisture content testing or in-situ monitoring data, preferably obtained by drying method or time domain reflectometer measurement. For areas without monitoring points, empirical estimation can be made based on topographic slope and historical rainfall data, and a continuous distribution field can be constructed using spatial interpolation methods (such as inverse distance weighting or kriging).

[0069] After completing the above parameter construction, consistency correction processing is performed on all parameter data. Specifically, the units of data from different sources are unified, the dimensions are normalized, and outliers are removed. Multi-source data at the same spatial location are fused to ensure the coordination between parameters. Subsequently, the processed parameter data is mapped to the numerical computing grid to establish a spatially discretized parameter field, where each grid cell corresponds to a complete set of unsaturated soil parameters.

[0070] In this embodiment, the computational unit refers to the basic grid unit formed after numerical discretization of the slope area, including internal computational units, surface boundary units, and potential slippage influence units. The boundary unit refers to the computational unit located at the slope surface, the exposed surface at the toe of the slope, the contact surface of the drainage ditch, or other boundaries where water exchange can occur, and is a specific type of computational unit. When calculating boundary flux, boundary water content, and boundary outward normal vector, b is used to represent the free drainage boundary unit index. When updating the internal permeability coefficient field, water content gradient response term, structural change response term, and seepage direction correction, i is used to represent the computational unit index.

[0071] The unsaturated parameter model data specifically includes: basic soil physical property parameters, water content and matrix suction characteristic curve data, unsaturated permeability coefficient function data, and initial stress and initial water content data;

[0072] The basic physical property parameters of the soil are used to characterize the soil structure and compaction.

[0073] The moisture content versus matrix suction characteristic curve data are used to characterize the unsaturated moisture retention properties;

[0074] In a preferred embodiment, the moisture content versus matrix suction characteristic curve data can be expressed using the Van Genuchten model, and the calculation formula is as follows:

[0075] ;

[0076] In the formula, It is the effective saturation. It is the volumetric moisture content. It is the residual moisture content. It is the saturated moisture content. It is matrix suction. , , These are the parameters for fitting the soil-water characteristic curve, preferably... This expression can be used to convert the results of indoor soil and water characteristic tests into functional parameters that can be directly called in subsequent seepage calculations.

[0077] After completing the soil-water characteristic curve fitting, the applicability of the fitting results is verified. Specifically, this includes checking whether the relationship between residual water content, saturated water content and measured water content meets the physical constraints; checking whether the fitted curve can cover the variation trend of the measured data points in the low suction region, medium suction region and high suction region respectively; and checking whether there are abrupt changes in the fitting parameters of adjacent soil layers that do not conform to the stratigraphic variation law.

[0078] For soil units with large fitting deviations, it is advisable to increase the number of test samples in the same layer or use the parameters of adjacent units for constraint correction. For soil units lacking complete test data, the initial parameter range can be selected according to the soil type, and the curve parameters can be corrected by subsequent field moisture content monitoring data.

[0079] The unsaturated permeability coefficient function data is used to characterize the seepage capacity under different water content conditions;

[0080] The initial stress and initial water-bearing state data are used as initial condition inputs for fluid-structure interaction analysis.

[0081] Example 3, see Figure 1 , Figure 2 This embodiment is based on the above embodiment. In step S2, the free drainage dynamic boundary evolution is used to characterize the dynamic change process of the boundary drainage capacity of the slope under free drainage conditions. Specifically, based on the water content of the surface area of ​​the slope, rainfall input information and water potential gradient distribution, the drainage capacity of the boundary area is dynamically determined, and a drainage response adjustment function related to time and local water content is constructed to continuously update the boundary drainage behavior and obtain dynamic boundary flux data.

[0082] The evolution of the free drainage dynamic boundary specifically adopts a free drainage dynamic boundary modeling method based on hysteresis response adjustment. In some preferred embodiments, a drainage response adjustment function can also be constructed to nonlinearly modulate the boundary drainage flux to reflect the hysteresis characteristics and gradual release characteristics in the drainage process.

[0083] The free drainage dynamic boundary modeling method based on hysteresis response adjustment includes the following steps:

[0084] Step S21: Boundary cell identification, used to determine the boundary regions participating in free drainage and establish the initial state. Specifically, the surface grid cells of the slope are screened to identify the boundary cells that meet the free drainage conditions, and the water content data, water potential data and rainfall input data corresponding to each boundary cell are extracted to obtain the initial boundary state data.

[0085] In a preferred embodiment, the boundary unit that meets the free drainage condition can be determined based on the slope's geometric exposure state, the direction of the water potential gradient, and the rainfall input state. Specifically, surface grid units located on the slope surface, the exposed surface at the toe of the slope, or the contact area with the drainage ditch are preferentially selected as candidate boundary units. It is then determined whether the water potential gradient of the candidate boundary unit has a component pointing outward from the boundary. When the candidate boundary unit is in the rainfall cessation stage, the weak rainfall stage, or the infiltration intensity is lower than the preset drainage initiation judgment value, and the volumetric water content of the boundary unit is higher than the critical water content for drainage initiation or there is an outward discharge water potential gradient in its neighborhood, the candidate boundary unit is identified as a free drainage boundary unit.

[0086] For boundary units that are still in the stage of heavy rainfall infiltration, they can be temporarily not treated as free drainage boundary units, but as rainfall infiltration boundary units. They can then enter the free drainage determination process after the rainfall input weakens or stops.

[0087] Step S22: Basic drainage flux calculation, used to obtain the basic drainage capacity of the boundary area under unsaturated conditions. Specifically, based on the unsaturated permeability coefficient function and the boundary water potential gradient, the basic drainage flux of each boundary unit is calculated and processed to obtain basic drainage flux data.

[0088] Step S23: Construction of drainage response regulation function, used to characterize the hysteresis and gradual release characteristics of drainage capacity during free drainage. Specifically, based on the water content data and time evolution information of the boundary unit, a nonlinear drainage response regulation function related to water content and time is constructed to modulate the drainage response process and obtain drainage response regulation function data.

[0089] The formula for calculating the drainage response adjustment function is as follows:

[0090] ;

[0091] In the formula, This is the drainage response adjustment function value of the boundary element at time t. As drainage response adjustment function data, it characterizes the degree of drainage release of the current boundary element under free drainage conditions. This is the volumetric water content of the current boundary cell, which can be directly obtained from the boundary state update result in step S2. It is the critical moisture content for drainage initiation, used to characterize the moisture content threshold for the transition of a boundary unit from a weakly drained state to a clearly drained state; It is the cumulative response time of the boundary unit since entering the current drainage stage, preferably from the moment the rainfall stops or from the moment the moisture content of the boundary unit exceeds the critical moisture content. The duration starting from the specified moment; It is a water content state sensitive adjustment coefficient, used to control the degree of influence of water content changes on the intensity of drainage response; It is the time response adjustment coefficient, used to control how quickly the drainage capacity is released over time;

[0092] Among them, the first item This is used to characterize the process by which the drainage capacity of a boundary unit gradually increases with changes in water content after the water content reaches the drainage initiation condition; when At this time, the boundary element can be considered to have not yet fully entered the effective drainage state, and the drainage response is relatively weak.

[0093] when Persistently greater than At that time, this item gradually increases with increasing moisture content and then tends to stabilize;

[0094] Second item This term is used to characterize the gradual release of drainage capacity over time, reflecting the lag and delay that actually exist in the free drainage process. In the early stage of drainage, this term increases slowly because the water redistribution and surface drainage channels have not yet been fully established. As time progresses, this term gradually increases, causing the drainage flux to transition from an initial suppressed state to a stable release state.

[0095] In some implementations, to avoid when When the exponential term exhibits unexpected amplification, (θ−θc) can be truncated to a non-negative value, i.e., when... At that time, take This ensures that the overall drainage response adjustment function remains within a reasonable range; furthermore, for ease of engineering implementation, the aforementioned The value of is preferably limited to between 0 and 1. When the calculation result exceeds the preset upper limit, it can be truncated to 1.

[0096] In the embodiments, if the slope soil is silty clay or a mixture of cohesive and silty soil, The optimal value is 0.18–0.30; for soils with high sand content and rapid drainage, The optimal value is 0.12–0.22; for soils with strong cohesion and high water retention capacity, The optimal value is 0.22 to 0.35. In the absence of specific test data, an initial value can be given based on the soil type, and then corrected by combining the drainage start time monitored on site.

[0097] In the embodiment, the A value of 2 to 8 can be preferred to control the sensitivity of moisture content changes to drainage initiation intensity; when the surface soil of the slope has a relatively uniform pore structure and the drainage initiation is relatively gentle, 2 to 4 can be preferred; when there is a significant loose layer on the surface of the slope and the drainage condition is sensitive to changes in water content, The preferred values ​​are 4 to 8, and in the preferred embodiment, the preferred values ​​are... ;

[0098] In the embodiment, the The preferred value is 0.05 to 0.50, with the unit varying accordingly. The time unit is selected and changes accordingly; when In hours The optimal value is 0.08–0.25; when there is a significant hysteresis effect in the boundary drainage, A value of 0.08–0.15 is acceptable; when drainage is established quickly, A value of 0.18 to 0.30 is acceptable; for typical scenarios of free drainage after rainfall, a value of [value missing] is preferred in the preferred embodiment. or ;

[0099] Furthermore, in an implementable example, when the current volumetric water content of a boundary cell is 0.28, the critical water content is... Take 0.20, Take 4.5, Pick And the cumulative response time after the unit enters the drainage stage. When the time is 6h, the corresponding drainage response adjustment function value can be calculated, which is used to nonlinearly modulate the basic drainage flux. Through this setting, the actual process of the boundary unit's drainage capacity being released slowly, then gradually increasing and stabilizing after the rainfall stops can be better reflected.

[0100] Step S24: Dynamic boundary flux correction, used to dynamically adjust the basic drainage flux. Specifically, based on the drainage response adjustment function data, the basic drainage flux is nonlinearly corrected, and the flux is smoothly updated in conjunction with time continuity constraints to obtain the drainage flux distribution data of the boundary unit. The specific calculation formula is as follows:

[0101] ;

[0102] ;

[0103] ;

[0104] In the formula, This is the initial corrected drainage flux of the b-th boundary cell at time t. The basic drainage flux is calculated based on the unsaturated permeability coefficient and water potential gradient. It is the value of the drainage response adjustment function. It is the actual boundary drainage flux after time continuity constraints. It is the actual boundary drainage flux of the b-th boundary element at the previous time t-1. This is the flux smoothing coefficient, preferably between 0.2 and 0.7. The basic unsaturated permeability coefficient is determined by the unsaturated permeability coefficient function in step S1. It is the volumetric water content of the boundary unit. It is the head gradient of the boundary unit. It is the outer normal vector of the boundary;

[0105] Step S25: Drainage capacity evolution update, used to construct the time change process of drainage capacity. Specifically, based on the time step advancement mechanism, the drainage flux of each boundary unit is continuously updated and processed to form the sequence data of drainage flux changing with time, and the drainage capacity changes with time data is obtained.

[0106] Step S26: Boundary water state response modeling, which is used to reflect the feedback effect of the drainage process on the boundary water state. Specifically, based on the change of drainage flux of the boundary unit, the water content of the boundary region is updated and calculated, and the trend of water state change over time is recorded to obtain the water state response data of the boundary region.

[0107] Preferably, the updated calculation of the water content of the boundary region can be based on the water balance of the boundary unit. For the b-th boundary unit, the current volumetric water content can be calculated based on the volumetric water content of the previous time step, the rainfall input in the current time step, the recharge of adjacent internal units, and the outward discharge flux. The calculation formula is as follows:

[0108] ;

[0109] In the formula, It is the volumetric water content of the b-th boundary unit at the next moment. It is the volume of the boundary element. It is the boundary drainage surface area. It is the amount of rainfall infiltration acting on the b-th boundary unit within the current time step. It represents the amount of water supplied from adjacent internal units to the b-th boundary unit within the current time step. It is the actual boundary drainage flux after time continuity constraints. It is the current time step. The total displacement is the outward displacement of the b-th boundary cell within the current time step. It is a function that limits the calculation results to the range of residual moisture content to saturated moisture content; when rainfall input or adjacent internal unit supply is obtained in the form of flux or flow, it can be converted into the water volume in the current time step according to the corresponding area of ​​action and time step and then substituted into the above formula;

[0110] Step S27: Update the hysteresis adjustment parameters to improve the adaptive capability of the drainage response model. Specifically, based on historical response data during the drainage process, the key parameters in the drainage response adjustment function are dynamically adjusted to obtain drainage hysteresis adjustment parameter data.

[0111] In a preferred embodiment, the historical response data includes a sequence of measured moisture content changes in the boundary region, drainage outflow monitoring data, the moment rainfall stopped, the start time of the boundary unit entering an effective drainage state, and a sequence of boundary flux changes in adjacent time steps.

[0112] When updating the lag adjustment parameters, priority is given to the time response adjustment coefficient. Corrections should be made; when the calculated peak drainage flux occurs significantly earlier than the measured peak drainage flux, the value should be reduced. The value of [value] is chosen to enhance drainage lag; when the calculated peak drainage flux significantly lags behind the measured peak drainage flux, [the value] is increased. The value of is chosen to accelerate the release of the drainage response;

[0113] Sensitive adjustment coefficient for water content The rate of change can be slightly adjusted based on the degree of matching between the rate of decrease in water content and the rate of increase in boundary flux; when the increase in drainage flux is insufficient after the water content reaches the critical water content, it can be appropriately increased. When the drainage flux is too sensitive to changes in water content and causes oscillations between adjacent time steps, the flow rate can be appropriately reduced. The and The updates are preferably limited to the aforementioned parameter value range, and the magnitude of a single update does not exceed 10% to 20% of the current parameter value;

[0114] The dynamic boundary flux data specifically includes boundary unit drainage flux distribution data, drainage capacity variation data over time, boundary region water content response data, and drainage hysteresis adjustment parameter data.

[0115] Example 4, see Figure 1 , Figure 3 This embodiment is based on the above embodiment. In step S3, the seepage path modeling is used to describe the evolution process of the seepage path inside the slope under the dynamic boundary drive. Specifically, based on the dynamic boundary flux data obtained in step S2, the water content distribution and its spatial gradient inside the slope are analyzed. Combined with the local deformation state or structural change information of the soil, the permeability coefficient field is adaptively updated, and the seepage path is directionally corrected and continuously reconstructed to obtain seepage path evolution data.

[0116] The permeability path modeling specifically introduces a path enhancement factor based on the water cut gradient to weight and adjust the local seepage direction, including the following steps:

[0117] Step S31: Internal state field extraction, used to obtain basic state information related to the evolution of seepage path inside the slope and determine potential seepage concentration areas. Specifically, based on the dynamic boundary flux data obtained in step S2, the water content distribution inside the slope is updated and calculated, and the spatial gradient analysis of water content is performed. At the same time, the local structural change index of each unit is extracted, and the areas that meet the water content gradient and structural change threshold conditions are screened to obtain seepage candidate area data.

[0118] In a preferred embodiment, the seepage candidate region data can be determined jointly based on the water content gradient response term and the structural change response term. For the i-th calculation unit, when its water content gradient response term... Greater than the preset gradient screening threshold, or structural change response item When the value exceeds the preset structural change screening threshold, the calculation unit is identified as a candidate unit; when neither of the two thresholds exceeds the corresponding threshold, but the weighted combination result of the two exceeds the comprehensive screening threshold, the calculation unit can also be identified as a candidate unit; the gradient screening threshold can preferably be 0.45 to 0.65, the structural change screening threshold can preferably be 0.30 to 0.50, and the comprehensive screening threshold can preferably be 0.50 to 0.70.

[0119] Step S32: Adaptive update of the permeability coefficient field, used to jointly adjust the permeability capacity according to the change in water content and the local structural response. Specifically, based on the unsaturated permeability coefficient function, the permeability coefficient of each unit is updated, and the water content gradient response term and the structural change response term are introduced to correct the permeability coefficient, so as to obtain the updated data of the spatial distribution of permeability coefficient.

[0120] The formula for calculating the spatial distribution update data of the permeability coefficient is:

[0121] ;

[0122] In the formula, It is the updated permeability coefficient of the i-th computational unit at time t, used to characterize the actual permeability of the unit under the current water content and local structural conditions; It is the i-th calculation unit at time t based on the current volumetric water content. The determined basic unsaturated permeability coefficient can be directly obtained from the unsaturated permeability coefficient function constructed in step S1 and used as the basic term; It is the water content gradient response weighting coefficient, used to control the influence of the water content gradient on the degree of permeability enhancement; It is the water content gradient response term of the i-th computational unit at time t, used to characterize the strength of water migration driving force in the region where the unit is located; It is the structural change response weighting coefficient, used to control the impact of local structural changes on the magnitude of permeability correction. It is the structural change response term of the i-th computational unit at time t, used to characterize the degree of change in permeability of the unit caused by local deformation, loosening or microcrack development;

[0123] Among them, basic items Used to reflect the original permeability of unsaturated soil under the current water content; First correction term This term reflects the more significant the moisture content gradient, the more pronounced the local water migration trend, thus enhancing permeability; the second correction term... This is used to reflect the potential expansion or reorganization of permeability channels under conditions of local structural relaxation, shear disturbance, or microcrack propagation, which can lead to enhanced permeability. For structural changes that are mainly characterized by compaction and closure, the response term of the structural change can be reduced or the structural enhancement correction can be avoided by misjudging the decrease in permeability caused by compaction as an increase in permeability.

[0124] By using the above multiplicative joint correction method, the combined influence of water migration drive and soil structure state on the permeability coefficient field can be considered simultaneously, making the evolution process of the seepage path more consistent with the actual change law under the action of free drainage boundary.

[0125] In a preferred embodiment, the moisture content gradient response term is calculated using a normalized gradient form, and the calculation formula is as follows:

[0126] ;

[0127] In the formula, It is the spatial gradient magnitude of water content in the i-th computational unit at time t. It is a preset gradient normalization upper limit, which can be determined based on historical operating conditions or the 95th percentile value of the gradient distribution across the entire field.

[0128] The structural change response term can be characterized by the local volumetric strain and shear strain of the previous time step, and the calculation formula is as follows:

[0129] ;

[0130] In the formula, It is the volumetric strain weighting coefficient. It is the local volumetric strain characteristic of the previous time step. It is a reference threshold for local volumetric strain characteristics. It is the shear strain weighting coefficient. It is a local shear strain characteristic. It is a reference threshold for local shear strain characteristics;

[0131] The local volumetric strain characteristics and local shear strain characteristics, for the initial time step, can be determined based on the initial stress field, geological fracture investigation results, or preset initial values ​​of structural disturbance.

[0132] and Preferred satisfaction When the main structural changes of the slope soil are shear disturbance, crack propagation, or local slippage, The preferred values ​​are 0.6 to 0.8. The optimal value is 0.2 to 0.4; when the slope soil exhibits volume relaxation, increased porosity, or unloading deformation as the main structural changes, The preferred values ​​are 0.5 to 0.7. The preferred value is 0.3 to 0.5;

[0133] and The value can be determined based on indoor triaxial tests, direct shear tests, or existing engineering experience. Alternatively, the 90%–95% quantile of the statistical distribution of volumetric strain and shear strain across the entire field at the corresponding moment can be used. For elements that are predominantly compressive and do not form conductive fractures, the effect of the volumetric strain term on the overall structure can be reduced. The contribution, or simply using shear strain terms and porosity increments to characterize structural changes;

[0134] In some embodiments, the moisture content gradient response term The calculation can be performed based on the difference in volumetric water content between the current cell and its neighboring cells, preferably using a normalized gradient form, ensuring its value is between 0 and 1; when When the moisture content is close to 0, it indicates that the moisture content distribution in the area is relatively uniform and the driving force of moisture migration is weak; when... When the value is close to 1, it indicates that there is a significant water content gradient in the area, which is prone to localized seepage concentration.

[0135] In some implementations, the structural change response term It can be obtained by normalizing local volumetric strain, shear strain, displacement increment, porosity change, or crack development index; preferably, it is expressed in the form of normalized absolute value of local volumetric strain; when When the value is small, it indicates that the local structural changes are not obvious; when When the value is large, it indicates that there has been a significant structural adjustment in the local soil, which may create more favorable conditions for seepage.

[0136] To ensure that the updated permeability coefficient is within a reasonable range, during project implementation, [the following can be done / performed]: and Set separate upper limits, for example, no greater than 2.5 and 2.0 respectively, or directly set the upper limit values. Set an overall upper limit so that it does not exceed the baseline permeability coefficient. The value is 3 to 6 times higher to avoid local outliers causing instability in subsequent numerical solutions.

[0137] In the embodiment, the A value of 0.5 to 3.0 can be preferred to control the enhancement of the permeability coefficient by the moisture content gradient; when the soil layers inside the slope are relatively homogeneous and local preferential flow is not obvious, The optimal value is 0.5 to 1.2; when there is a significant wetting front advancing or local gradient accumulation phenomenon inside the slope, The preferred value is 1.2 to 2.5; in the preferred embodiment, it can be 1.2 to 2.5. or ;

[0138] In the embodiment, the A value of 0.3 to 2.0 can be preferred to control the degree of correction of permeability by structural changes; when the soil structure changes are small or the slope is in the early response stage, The preferred value is 0.3 to 0.8; when local deformation is significant, microcracks are propagating, or the structure tends to loosen, The preferred value is 0.8 to 1.8; in the preferred embodiment, it can be taken as... or ;

[0139] Furthermore, in an implementable example, if the current base unsaturated permeability coefficient of a certain computing unit... for Moisture content gradient response term The structural change response term is 0.60. It is 0.35. Take 1.5, If we take 1.0, the updated permeability coefficient can be directly calculated using the above formula; this result can be used to characterize the enhanced permeability effect of the unit at the current moment due to the large local water content gradient and certain structural changes.

[0140] Step S33: Seepage direction correction, used to enhance the guiding effect of the water content gradient on the seepage path direction. Specifically, a path enhancement factor is constructed based on the water content gradient of each unit, and the path enhancement factor is weighted and fused with the original seepage direction to correct the seepage direction, thereby obtaining seepage path direction distribution data. The calculation formula is as follows:

[0141] ;

[0142] ;

[0143] In the formula, This is the corrected dominant seepage direction, i.e., the distribution data of seepage path direction. It is path enhancement weight. The original seepage direction is determined by the head gradient. It is the path enhancement direction determined by the moisture content gradient. This represents the maximum enhancement weight, preferably between 0.6 and 0.8. This is the gradient response scaling factor;

[0144] Step S34: Local priority flow channel identification, used to identify continuous channel structures where seepage is concentrated within the slope. Specifically, based on the updated spatial distribution data of permeability coefficient and the distribution data of seepage path direction, the continuity of seepage between adjacent units is determined, and units that meet the conditions of enhanced permeability and directional consistency are connected and merged to obtain local priority flow channel identification data.

[0145] In a preferred embodiment, the seepage continuity between adjacent units can be determined based on the degree of permeability enhancement and the consistency of seepage direction. For adjacent i-th and j-th calculation units, when both belong to the seepage candidate region, and the enhancement factor of their updated permeability coefficient relative to the basic unsaturated permeability coefficient is not lower than the preset enhancement threshold, and the angle between their corrected seepage dominant directions is less than the preset direction angle threshold, it is determined that there is seepage path continuity between the two calculation units.

[0146] The enhancement threshold can preferably be set to 1.3–2.0, and the directional angle threshold can preferably be set to 30°–45°. For slopes with obvious fissure flow characteristics, the directional angle threshold can be appropriately relaxed; for homogeneous soil slopes or conditions with high monitoring noise, the enhancement threshold can be appropriately increased. Subsequently, adjacency search or connected component marking methods can be used to connect and merge adjacent units that meet the continuity determination conditions to form one or more local priority flow channels. To avoid misjudgment caused by isolated abnormal units, channel segments with fewer than a preset number of consecutive units can be removed. The preset number can preferably be set to 3–5 adjacent units.

[0147] Step S35: Seepage concentration area identification, used to determine seepage concentration areas and construct the continuous evolution results of seepage paths. Specifically, based on the local priority flow channel identification data, the degree of seepage concentration in the area surrounding the channel is analyzed, and areas that meet the concentration conditions are identified. At the same time, the seepage paths are continuously reconstructed by combining time series information to obtain seepage concentration area identification data.

[0148] In a preferred embodiment, the area surrounding the channel can be determined based on the adjacency relationship of each unit in the local priority flow channel. Preferably, the first-order or second-order neighboring units adjacent to the channel unit are taken as the evaluation area. For the calculation unit in the evaluation area, its permeability enhancement factor, water content gradient response term, boundary flux transfer direction and spatial distance from the identified channel can be considered to determine whether it belongs to the seepage concentration area.

[0149] When a unit is adjacent to a local preferential flow channel and its water content gradient response term or permeability enhancement factor is consistently higher than the corresponding threshold, it is identified as a seepage concentration unit; when this state remains stable for more than two consecutive time steps, it is identified as a persistent seepage concentration region.

[0150] Isolated high-value units that appear only in a single time step can be recorded as transient anomaly units, but are not directly output as stable seepage concentration areas.

[0151] Step S36: Output the permeability path evolution results, which is used to structure and organize the results of each sub-step and form a unified output. Specifically, it involves associating and integrating the spatial distribution update data of the permeability coefficient, the distribution data of the permeation path direction, the identification data of local preferential flow channels, and the identification data of the permeation concentration area to obtain the permeability path evolution data.

[0152] When outputting the seepage path evolution data, it is preferable to use the calculation unit as the basic index and record the updated permeability coefficient, dominant seepage direction, water content gradient response term, structural change response term, whether it belongs to a local preferential flow channel, whether it belongs to a seepage concentration area, and the corresponding time step information for each unit. For units identified as local preferential flow channels, their channel number, channel length, number of adjacent connected units, and channel duration can also be recorded.

[0153] The seepage path evolution data specifically includes: permeability coefficient spatial distribution update data, seepage path direction distribution data, local preferential flow channel identification data, and seepage concentration area identification data;

[0154] The spatial distribution update data of the permeability coefficient is used to characterize the permeability capacity distribution of each unit inside the slope at the current moment;

[0155] The seepage path direction distribution data is used to characterize the dominant seepage direction of each unit after direction correction;

[0156] The local preferential flow channel identification data is used to characterize the high-permeability continuous channel structure formed inside the slope;

[0157] The data identifying concentrated seepage areas is used to characterize key areas where local seepage accumulation is significant and may induce subsequent response changes.

[0158] By performing the above operations, this solution addresses the technical problem in existing unsaturated seepage analysis methods that have a fixed permeability coefficient field or rely solely on the single variable of water content, making it difficult to reflect local structural changes and the formation process of preferential flow channels. This solution creatively adopts an adaptive update model of permeability coefficient based on the response of water content gradient and structural change, and combines it with a path enhancement factor to correct the seepage direction, while simultaneously identifying and continuously reconstructing local preferential flow channels.

[0159] For example, when local shear deformation or microcrack development occurs on a slope, traditional models still use the original permeability coefficient, which makes it impossible to capture the evolution process of local seepage concentration and channel formation leading to rapid infiltration. However, this scheme introduces water content gradient response terms and structural change response terms to multiplicatively correct the permeability, so that the permeability coefficient increases synchronously with water content gradient and structural changes. This enables the identification of key risk areas such as concentrated seepage zones at the toe of the slope or fractured flow paths.

[0160] Example 5, see Figure 1 This embodiment is based on the above embodiment. In step S4, the fluid-structure interaction response modeling is used to perform coupled modeling of the mechanical response of the slope soil based on the seepage field update. Specifically, based on the seepage path evolution data obtained in step S3, the distribution of pore water pressure inside the soil is updated, and the stress field is calculated and processed in combination with the strength characteristics of unsaturated soil to realize the linkage update between the seepage field and the stress field and obtain fluid-structure interaction response data.

[0161] In this embodiment, the fluid-structure interaction response modeling can be implemented in the following way:

[0162] First, based on the seepage path evolution data obtained in step S3, the updated spatial distribution data of the permeability coefficient and the distribution data of the seepage path direction are mapped to the computational grid as input parameters for the current seepage field calculation. On this basis, the pore water pressure inside the slope is solved and updated by combining the mass conservation equation and Darcy's law to obtain the pore water pressure distribution of each computational unit at the current moment.

[0163] Subsequently, based on the updated pore water pressure distribution data and combined with the effective stress principle, the soil stress field was corrected. The effective stress of the soil can be obtained by subtracting the pore water pressure from the total stress. On this basis, a conventional unsaturated soil constitutive model, such as the Mohr-Coulomb model or an elastoplastic constitutive model, was selected to solve the stress-strain relationship of each calculation unit, and the soil stress and strain response data were obtained.

[0164] Furthermore, based on the strain results obtained from the calculation, the displacement field of the overall slope and local areas is updated and calculated. The finite element method or finite difference method is preferably used to iteratively solve the nodal displacements, thereby obtaining local deformation and displacement field data.

[0165] Regarding the updating of strength parameters, the cohesion parameters or shear strength parameters can be modified based on the expression relationship of unsaturated soil strength. For example, the strength parameters can be adjusted in segments or updated by function mapping according to the current matrix suction or water content state, so that they can reflect the changes in the mechanical properties of the soil under different water content states.

[0166] The fluid-structure interaction response data specifically includes: pore water pressure distribution data, soil stress and strain response data, local deformation and displacement field data, and strength parameter update data.

[0167] Example 6, see Figure 1 This embodiment is based on the above embodiment. In step S5, the fluid-structure interaction solution control is used to perform numerical solution and calculation control on the fluid-structure interaction process. Specifically, based on the fluid-structure interaction response data obtained in step S4, the time step, iteration convergence condition and error control strategy in the calculation process are adjusted and processed, and multi-step iterative solution is performed to obtain slope stability analysis result data.

[0168] The fluid-structure interaction solution control can be implemented in this embodiment as follows:

[0169] First, based on the fluid-structure interaction response data obtained in step S4, a time step advancement mechanism is established to control the fluid-structure interaction process in steps. In each time step, the initial calculation conditions are first set according to the current seepage field and stress field state, and the seepage equation and mechanical equation are coupled and solved according to the preset time step.

[0170] Secondly, in terms of time step control, the time step can be adaptively adjusted according to the change in pore water pressure or displacement increment. When the change in calculation results is relatively gentle, the time step can be appropriately increased, and when the change is more drastic, the time step can be reduced to improve the stability and efficiency of the calculation.

[0171] In terms of iterative convergence control, pore water pressure residual, displacement residual, and stress balance error can be set as convergence criteria. When all errors are less than a preset threshold, the current time step is considered converged; otherwise, iterative updates continue. Preferably, the residual threshold can be set to 10. −3 ~10 −5 Magnitude;

[0172] In terms of error control strategies, the calculation results between consecutive time steps can be smoothed to avoid the adverse effects of local numerical oscillations on the overall calculation results. At the same time, abnormal units, such as abrupt changes in permeability coefficient or strain, can be marked and corrected in subsequent iterations.

[0173] After completing all time step calculations, the fluid-structure interaction results at each time point are summarized and analyzed. Based on the soil strength parameters and stress distribution, the slope stability is evaluated. Preferably, the slope safety factor can be calculated using the strength reduction method or the limit equilibrium method, and potential slip surfaces can be identified by combining the displacement field and strain concentration areas.

[0174] The slope stability analysis results data specifically include: slope stability evaluation index data, safety factor calculation results data, potential slip surface distribution data, and instability risk area identification data.

[0175] By performing the above operations, this solution addresses the technical problems of existing fluid-structure interaction numerical analysis methods, such as lag in the coupling update of seepage field and stress field, insufficient stability of the solution process, and weak engineering applicability. It creatively adopts a fluid-structure interaction response modeling and adaptive solution control mechanism driven by seepage path evolution. By dynamically updating pore water pressure, stress, strain, and strength parameters, and introducing adaptive time step adjustment, residual convergence control, and abnormal element correction strategies, it achieves stable and efficient solution of the fluid-structure interaction process.

[0176] For example, during rainfall-induced slope deformation, rapid changes in pore pressure often lead to numerical divergence or calculation oscillation. This scheme dynamically adjusts the time step based on the magnitude of pore pressure changes and sets multiple residual criteria. At the same time, it marks and corrects abrupt permeability changes, thereby ensuring stable convergence even in the strongly nonlinear stage and improving the reliability and applicability of engineering calculations.

[0177] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

[0178] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention.

[0179] The present invention and its embodiments have been described above. This description is not restrictive, and the accompanying drawings are only one embodiment of the present invention; the actual structure is not limited thereto. In conclusion, if those skilled in the art are inspired by this description and design similar structures and embodiments without departing from the spirit of the invention, such designs should fall within the protection scope of the present invention.

Claims

1. A fluid-structure interaction numerical analysis method for unsaturated slopes under free drainage conditions, characterized in that: The method includes the following steps: Step S1: Unsaturated multi-source parameter modeling to obtain unsaturated parameter model data; Step S2: Dynamic boundary evolution of free drainage. Based on the water content of the slope surface area, rainfall input information and water potential gradient distribution, the drainage capacity of the boundary area is dynamically determined, and a drainage response adjustment function related to time and local water content is constructed to continuously update the boundary drainage behavior and obtain dynamic boundary flux data. Step S3: Seepage path modeling. Based on dynamic boundary flux data, the distribution of water content and spatial gradient inside the slope are analyzed. Combined with the local deformation state or structural change information of the soil, the permeability coefficient field is adaptively updated, and the seepage path is directionally corrected and continuously reconstructed to obtain seepage path evolution data. In step S3, the permeation path modeling specifically introduces a path enhancement factor based on the water content gradient to weight and adjust the local seepage direction, including the following steps: internal state field extraction, adaptive update of the permeability coefficient field, correction of seepage direction, identification of local preferential flow channels, identification of seepage concentration areas, and output of permeation path evolution results. The internal state field extraction is based on the dynamic boundary flux data obtained in step S2. The internal water content distribution of the slope is updated and calculated, and the spatial gradient of water content is analyzed. At the same time, the local structural change index of each unit is extracted, and the areas that meet the water content gradient and structural change threshold conditions are screened to obtain seepage candidate area data. The adaptive update of the permeability coefficient field is based on the unsaturated permeability coefficient function. The permeability coefficient of each unit is updated, and the water content gradient response term and the structural change response term are introduced to correct the permeability coefficient, so as to obtain the permeability coefficient spatial distribution update data. The seepage direction correction is achieved by constructing a path enhancement factor based on the water content gradient of each unit, and then weighting and fusing the path enhancement factor with the original seepage direction to correct the seepage direction and obtain seepage path direction distribution data. The local priority flow channel identification is based on the updated spatial distribution data of permeability coefficient and the distribution data of seepage path direction. It determines the seepage continuity between adjacent units and connects and merges units that meet the conditions of enhanced permeability and directional consistency to obtain local priority flow channel identification data. The seepage concentration area identification is based on the local priority flow channel identification data. The degree of seepage concentration in the area surrounding the channel is analyzed, and the areas that meet the concentration conditions are identified. At the same time, the seepage path is continuously reconstructed by combining time series information to obtain seepage concentration area identification data. The output of the permeation path evolution results involves associating and integrating the permeability coefficient spatial distribution update data, permeation path direction distribution data, local preferential flow channel identification data, and permeation concentration area identification data to obtain permeation path evolution data. Step S4: Fluid-structure interaction response modeling. Based on the seepage path evolution data, the distribution of pore water pressure inside the soil is updated, and the stress field is calculated and processed in combination with the strength characteristics of unsaturated soil to obtain fluid-structure interaction response data. Step S5: Fluid-structure interaction solution control. Based on the fluid-structure interaction response data, the time step, iteration convergence conditions, and error control strategies in the calculation process are adjusted to obtain the slope stability analysis results data.

2. The numerical analysis method for fluid-structure interaction under free drainage conditions of an unsaturated slope according to claim 1, characterized in that: In step S1, the unsaturated parameter model data specifically includes: basic soil physical property parameters, water content and matrix suction characteristic curve data, unsaturated permeability coefficient function data, and initial stress and initial water content data.

3. The fluid-structure interaction numerical analysis method for unsaturated slopes under free drainage conditions according to claim 2, characterized in that: In step S2, the evolution of the free drainage dynamic boundary specifically adopts a free drainage dynamic boundary modeling method based on hysteresis response adjustment, including the following steps: boundary unit identification, basic drainage flux calculation, drainage response adjustment function construction, dynamic boundary flux correction, drainage capacity evolution update, boundary water state response modeling, and hysteresis adjustment parameter update. The boundary unit identification process involves filtering the surface grid units of the slope to identify boundary units that meet the free drainage conditions, and extracting the water content data, water potential data, and rainfall input data corresponding to each boundary unit to obtain the initial boundary state data. The basic drainage flux calculation is based on the unsaturated permeability coefficient function and the boundary water potential gradient. The basic drainage flux of each boundary unit is calculated and processed to obtain basic drainage flux data. The drainage response regulation function is constructed based on the water content data and time evolution information of the boundary cells. A nonlinear drainage response regulation function related to the water content state and time is constructed to modulate the drainage response process and obtain drainage response regulation function data.

4. The fluid-structure interaction numerical analysis method for unsaturated slopes under free drainage conditions according to claim 3, characterized in that: In step S2, the dynamic boundary flux correction is performed by nonlinearly correcting the basic drainage flux based on the drainage response adjustment function data, and by smoothly updating the flux in combination with time continuity constraints, to obtain the drainage flux distribution data of the boundary unit. The drainage capacity evolution update is based on a time-step advancement mechanism, which continuously updates the drainage flux of each boundary unit and forms a sequence of drainage flux changing over time, thus obtaining drainage capacity changing over time data. The boundary water state response modeling is based on the change in drainage flux of the boundary unit, updating the water content of the boundary region, and recording the trend of water state change over time to obtain the water state response data of the boundary region. The hysteresis adjustment parameter update is based on historical response data during the drainage process, and the key parameters in the drainage response adjustment function are dynamically adjusted to obtain drainage hysteresis adjustment parameter data.

5. The fluid-structure interaction numerical analysis method for unsaturated slopes under free drainage conditions according to claim 4, characterized in that: In step S4, the fluid-structure interaction response data specifically includes: pore water pressure distribution data, soil stress and strain response data, local deformation and displacement field data, and strength parameter update data.

6. The fluid-structure interaction numerical analysis method for unsaturated slopes under free drainage conditions according to claim 5, characterized in that: In step S5, the slope stability analysis results data specifically include: slope stability evaluation index data, safety factor calculation results data, potential slip surface distribution data, and instability risk area identification data.

Images