Stratum rainfall seepage deformation coupling numerical simulation method, system and storage medium

Abstract

The present application relates to geotechnical engineering disaster monitoring and early warning technical field, disclose stratum rainfall seepage deformation coupling numerical simulation method, system, storage medium, the method comprises: the three-dimensional geological model reflecting the characteristics of soft and hard rock interbed is established;Real-time conversion of field rainfall data into hydraulic boundary conditions;Calculate the shear dilatancy of rock mass and determine the chemical damage acceleration factor, update the total damage variable;According to the spatial gradient of damage variable, the permeability tensor representing the dominant flow channel is reconstructed;Based on the permeability tensor, the fluid-solid coupling equation is solved, and the pore water pressure and the calculated displacement field are obtained;Using the field monitoring data to construct the likelihood function to inverse modify the model parameters;Based on the corrected calculation results, the stability is evaluated and the early warning signal is output. The present application realizes the accurate simulation of the landslide evolution process under complex geological environment by constructing the force coupling mechanism and the dominant flow channel model.

Description

Technical Field

[0001] This invention relates to the field of geotechnical engineering disaster monitoring and early warning technology, specifically to a method, system, and storage medium for coupled numerical simulation of ground rainfall, seepage, and deformation. Background Technology

[0002] The Wumeng Mountain region has a complex geological structure, with the Feixianguan Formation widely distributed throughout the area. This stratum exhibits a typical interbedded structure of soft and hard rocks, with alternating layers of mudstone and hard rock. Due to the region's humid climate and frequent, intense rainfall, the rock mass on slopes is prone to alteration due to rainwater infiltration. This unique geological structure, combined with hydrological and meteorological conditions, leads to frequent landslides and other geological disasters in the region, posing a severe challenge to the safety of local engineering construction and operation.

[0003] For the analysis of such geological hazards, current engineering practices mostly employ numerical simulation techniques for assessment. The common approach involves constructing a generalized geological model based on geological survey data, simplifying rainfall into flow rate or hydraulic head boundaries applied to the model surface. By calculating the changes in pore water pressure during rainfall infiltration, and combining this with soil and rock strength criteria, the stress field distribution is analyzed. Finally, a safety factor is calculated based on limit equilibrium theory or the strength reduction method to determine the slope's stability under rainfall conditions.

[0004] However, existing technologies have shortcomings in handling such complex coupled problems. First, traditional methods are mostly based on the assumption of continuum mechanics, making it difficult to accurately describe discontinuous large deformation behaviors such as shear slippage and normal opening at the interface between soft and hard rocks in the Feixianguan Formation. Second, regarding multi-field coupling mechanisms, the nonlinear accelerating effect of mechanical shear expansion on chemical dissolution is often ignored, failing to reflect the guiding mechanism of damage gradient on the formation of dominant flow channels, resulting in a disconnect between the simulation of high-permeability channels and actual fracture evolution. In addition, model parameters are mostly taken from static laboratory tests, lacking a mechanism for dynamic inversion and correction of deep mechanical coupling parameters using real-time displacement data from the field, making it difficult for the calculation results to truly reflect the time-varying damage state at deep formations. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a numerical simulation method, system, and storage medium for coupled formation rainfall, seepage, and deformation. This solves the problem that existing technologies cannot accurately simulate the complex instability process of the Feixianguan Formation in the Wumeng Mountains, where interbedded soft and hard rocks are combined with interlayer slippage, shear expansion, accelerated chemical damage, and the evolution of dominant flow channels under rainfall conditions, resulting in insufficient accuracy in landslide disaster prediction.

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] On the one hand, this invention provides a numerical simulation method for coupled formation rainfall, seepage, and deformation, comprising the following steps:

[0008] A three-dimensional geological model with interlayered soft and hard rocks was established, and the rainfall data collected on site was converted into the hydraulic boundary conditions of the three-dimensional geological model in real time.

[0009] Calculate the shear dilatation rate of rock mass units in the three-dimensional geological model, thereby determining the chemical damage acceleration factor, and then update the total damage variable of the rock mass units through the chemical damage acceleration factor.

[0010] The permeation tensor characterizing the dominant flow channel is reconstructed based on the distribution gradient of the total damage variable in three-dimensional space; then, the fluid-structure interaction equation is solved based on the permeation tensor to obtain the pore water pressure distribution, thereby updating the effective stress and obtaining the calculated displacement field.

[0011] A likelihood function is constructed based on the deviation between the calculated displacement field and the displacement monitoring data collected on site. The model parameters of the three-dimensional geological model are then corrected by inversion using the likelihood function.

[0012] The modified three-dimensional geological model outputs the rainfall-seepage-deformation coupling results, thereby assessing formation stability and outputting early warning signals.

[0013] By employing a coupling mechanism that correlates the volumetric expansion effect caused by rock mass shear deformation with the rate of chemical dissolution damage, and introducing a permeability tensor reconstruction technique based on damage gradient, this method achieves a realistic reflection of the multi-field coupled evolution of force, water, and chemical processes in the Feixianguan Formation of the Wumeng Mountains under rainfall conditions. This method overcomes the shortcomings of traditional methods that neglect the accelerating effect of rock mass shear fracture on chemical dissolution and the distortion in simulating dominant flow channels, thus improving the physical realism and spatiotemporal accuracy of landslide disaster prediction under complex geological conditions.

[0014] Preferably, establishing a three-dimensional geological model reflecting the interbedded soft and hard rocks of the Feixianguan Formation in the Wumeng Mountains includes: processing borehole exploration data using an interpolation algorithm to construct stratigraphic interfaces; performing geometric cutting operations on the stratigraphic interfaces to generate mudstone soft layer entities and hard rock layer entities; discretizing the mudstone soft layer entities and the hard rock layer entities to generate rock mass units in the three-dimensional geological model; identifying the contact boundaries between the mudstone soft layer entities and the hard rock layer entities; setting interface units with independent constitutive properties at the contact boundaries; setting the mechanical parameters of the interface units; and simulating the shear displacement and slip deformation behavior of the strata along the bedding planes through the interface units.

[0015] By introducing interface elements between soft and hard rock entities, the unique interlayer slippage and separation behavior of the Feixianguan Formation can be accurately simulated, solving the problem that continuous medium models cannot describe large deformations of discontinuous layers, thus more accurately capturing landslide failure modes initiated along weak interlayers.

[0016] Preferably, the process of converting the on-site collected rainfall data into hydraulic boundary conditions for a three-dimensional geological model in real time includes: identifying the upper surface boundary nodes of the three-dimensional geological model; extracting the real-time rainfall intensity from the on-site rainfall data and comparing it with the infiltration capacity of the surface soil at the upper surface boundary nodes; when the real-time rainfall intensity does not exceed the infiltration capacity, directly applying the real-time rainfall intensity as a second-type flow boundary condition; when the real-time rainfall intensity exceeds the infiltration capacity, modifying the hydraulic boundary conditions to a first-type constant head boundary condition, and simulating surface water infiltration through the first-type constant head boundary condition.

[0017] By dynamically determining the relationship between rainfall intensity and infiltration capacity and switching boundary condition types, the process of excess infiltration runoff and surface water infiltration under heavy rainfall conditions is realistically simulated. This avoids the situation where a single flow boundary leads to model calculation divergence or physical distortion under extreme rainfall conditions, thus improving the accuracy of hydraulic boundary conditions.

[0018] Preferably, calculating the shear expansion coefficient includes: using an elastoplastic constitutive algorithm to determine whether the rock mass element has entered the plastic yield state; for the rock mass element that has entered the plastic yield state, quantifying the volume expansion component caused by shear slip according to the non-associated flow rule; extracting the plastic volumetric strain rate corresponding to the volume expansion component, and defining the plastic volumetric strain rate as the shear expansion coefficient used to characterize the shear expansion strength.

[0019] The shear expansion behavior of rock mass after yielding was quantified using the non-associated flow law, and the plastic volumetric strain rate was accurately extracted as a physical index, providing a quantitative basis consistent with rock mechanics mechanism for establishing the coupling relationship between mechanical deformation and accelerated chemical damage.

[0020] Preferably, updating the total damage variable of the rock mass unit includes: constructing a calculation logic for a chemical damage acceleration factor, defining the chemical damage acceleration factor as an exponential function of a natural constant, wherein the exponent of the exponential function is the product of the mechanical coupling coefficient and the shear expansion rate; calculating the chemical damage evolution rate, wherein the chemical damage evolution rate is the product of a baseline chemical dissolution rate constant, a term reflecting the remaining integrity of the rock mass, and the chemical damage acceleration factor; calculating the chemical damage increment based on the chemical damage evolution rate, and updating the total damage variable by superimposing the chemical damage increment with the mechanical damage increment.

[0021] By constructing an exponential acceleration factor, the physical process of rock mass shear microcrack propagation leading to an increase in reaction surface area and subsequently an exponentially accelerated chemical dissolution is mathematically described. This mechanistically accelerated unidirectional strong coupling mechanism not only quantifies the nonlinear superposition effect of mechanical and chemical damage but also reveals the microscopic dynamic root of rock mass strength attenuation during rainfall-induced landslides.

[0022] Preferably, the reconstructing of the permeation tensor representing the dominant flow channel includes: calculating the gradient vector of the total damage variable at the center point of the grid cell in three-dimensional space using a numerical differential algorithm; constructing a local orthogonal coordinate system based on the direction of the gradient vector, setting the tangential permeation coefficient to grow nonlinearly exponentially with the total damage variable in the local orthogonal coordinate system; constructing a coordinate rotation matrix, and using the coordinate rotation matrix to transform the permeation attribute in the local orthogonal coordinate system back to the global coordinate system to obtain the permeation tensor representing the dominant flow channel in full tensor form.

[0023] An anisotropic permeability tensor is constructed using the spatial gradient of the damage field, realizing the mapping from scalar damage to tensor permeability. This method can automatically identify and simulate dominant flow channels (such as fracture networks) that dynamically expand with damage evolution, enabling the groundwater flow field to adaptively concentrate in high-damage areas, thus realistically reflecting the non-uniform seepage characteristics in fractured rock masses.

[0024] Preferably, solving the fluid-structure interaction equation based on the permeation tensor includes: calculating the local Reynolds number of the fluid by combining the characteristic length of the mesh element and the fluid viscosity, comparing the local Reynolds number with a preset critical Reynolds number threshold; when the local Reynolds number is lower than the critical Reynolds number threshold, solving the fluid-structure interaction equation using the linear Darcy's law; when the local Reynolds number is higher than or equal to the critical Reynolds number threshold, switching to the Forchheimer nonlinear permeation algorithm to solve the equation, and introducing an inertial drag term to correct the pressure gradient in the fluid-structure interaction equation.

[0025] By introducing the Reynolds number criterion to automatically switch the seepage control equations, the problem that high-velocity nonlinear seepage (turbulence) in the dominant flow channel cannot be described by Darcy's law is effectively solved. The introduction of an inertial drag term for correction improves the accuracy of pore water pressure distribution calculations under fracture inrush or high-flow-rate seepage conditions.

[0026] Preferably, constructing and correcting the model parameters of the three-dimensional geological model through the likelihood function includes: obtaining measured values ​​of surface displacement and deep displacement to form an observation vector; extracting the values ​​of the calculated displacement field at the corresponding locations to form a simulation vector; assuming that the observation error follows a multidimensional independent Gaussian distribution, constructing the likelihood function to quantify the degree of matching based on the residuals of the observation vector and the simulation vector; generating candidate samples based on the likelihood function using a Markov chain Monte Carlo sampling algorithm; and determining the corrected model parameters through statistical analysis. The model parameters include a baseline chemical dissolution rate constant and a mechanochemical coupling coefficient.

[0027] By using the Markov chain Monte Carlo method under the Bayesian inference framework, multi-source displacement monitoring data is integrated to invert and correct deep mechanical coupling parameters that are difficult to measure directly. This reduces the uncertainty of model parameters, making the calculation results gradually approach the field measured data over time, and improving the model's consistency with reality.

[0028] Preferably, the assessment of formation stability and output of early warning signals includes: identifying potential interlayer slip surface paths; calculating the proportion of the unit length exceeding the critical failure value on the interlayer slip surface path to the total length of the slip surface as the penetration damage degree; acquiring the displacement rate of key monitoring points and performing normalization processing; weighting and fusing the penetration damage degree and the normalized displacement rate to generate a landslide initiation risk index; and generating the early warning signal containing the potential slip path and the expected failure time when the landslide initiation risk index exceeds a preset threshold.

[0029] By integrating the through-damage state of structural surfaces with kinematic displacement rate information, a comprehensive landslide initiation risk index was constructed. Compared with a single displacement threshold or safety factor, this criterion can more comprehensively characterize the evolution process of landslides from local damage to overall instability, thus providing more forward-looking and reliable disaster early warning.

[0030] On the other hand, the present invention also provides a formation rainfall-seepage-deformation coupled numerical simulation system for performing the formation rainfall-seepage-deformation coupled numerical simulation method described in the foregoing aspect, comprising:

[0031] The data acquisition and monitoring module is used to obtain physical quantity inputs of the geological environment, and to collect static geological survey data and dynamic real-time monitoring data.

[0032] The geological modeling module is used to construct a three-dimensional geological grid model containing interface units based on the data input from the data acquisition and monitoring module.

[0033] The constitutive coupling calculation module is used to calculate the shear dilatancy rate and chemical damage acceleration factor of rock mass elements, and update the total damage variable.

[0034] The permeation channel reconstruction module is used to reconstruct the permeation tensor characterizing the dominant flow channel based on the distribution gradient of the total damage variable in three-dimensional space.

[0035] The fluid-structure interaction solution module is used to solve the fluid-structure interaction equation based on the permeation tensor characterizing the dominant flow channel to obtain the pore water pressure distribution;

[0036] The parameter inversion correction module is used to construct a likelihood function using field displacement monitoring data and perform inversion correction on the parameters in the constitutive coupling calculation module.

[0037] The early warning visualization module is used to receive calculation results and generate landslide disaster early warning signals based on stability criteria.

[0038] The aforementioned system integrates complex force-coupled algorithms and numerical calculation processes into a modular system, realizing automated full-process processing from data acquisition, model updating, coupled calculation to early warning release, providing an efficient computing platform and tool support for real-time monitoring and emergency decision-making of geological disasters in the Wumeng Mountain area.

[0039] Accordingly, the present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the described numerical simulation method for coupling formation rainfall, seepage, and deformation.

[0040] In summary, the present invention provides one or more technical solutions, which, compared with the prior art, have the following beneficial effects:

[0041] 1. This invention establishes interface units with independent constitutive properties at the interface between soft mudstone layers and hard rock layers, allowing for interlayer shear displacement and normal opening. This technical feature effectively solves the problem that traditional continuous medium models struggle to describe the significant discontinuous large deformation behavior in the interbedded soft and hard rock strata of the Feixianguan Formation in the Wumeng Mountains. It can accurately simulate the shear slippage and separation mechanisms of strata along weak interlayers, thereby improving the physical realism of landslide failure mode analysis for this specific geological structure.

[0042] 2. This invention determines the chemical damage acceleration factor by calculating the shear expansion rate of rock mass units and reconstructs the permeability tensor characterizing the dominant flow channel based on the spatial gradient of the total damage variable. These technical features establish a strong coupling mechanism of mechanically induced chemical corrosion, chemically eroded water, and weak water forces, quantifying the accelerating effect of mechanical crack propagation on chemical dissolution and the guiding effect of damage on the seepage path. This overcomes the shortcomings of simplified or disconnected multi-field coupling mechanisms in existing technologies, and improves the simulation accuracy of long-term strength decay and non-uniform seepage evolution processes in rock masses under rainfall conditions.

[0043] 3. This invention constructs a likelihood function by collecting on-site displacement monitoring data and uses the Markov chain Monte Carlo algorithm to invert and correct key model parameters such as the mechanical coupling coefficient. This technical feature feeds real-time monitoring data back into the numerical model, achieving dynamic calibration from theoretical model to on-site measurement. This effectively reduces the uncertainty caused by the difficulty in obtaining deep geological parameters, making the output stratum stability assessment results and landslide disaster early warning signals more closely match the actual on-site conditions, thus enhancing the timeliness and reliability of the early warning. Attached Figure Description

[0044] Figure 1 This is a schematic diagram of the overall architecture of the formation rainfall-seepage-deformation coupled numerical simulation system according to an embodiment of the present invention;

[0045] Figure 2 This is a flowchart illustrating the numerical simulation method for coupled formation rainfall, seepage, and deformation according to an embodiment of the present invention.

[0046] Figure 3 This is a time history curve of cumulative displacement and interlayer damage evolution at key monitoring points of a slope during a heavy rainfall event, as output by the system in this embodiment of the invention.

[0047] The module includes: 10. Data acquisition and monitoring module; 20. Geological modeling module; 30. Constitutive coupling calculation module; 40. Permeability channel reconstruction module; 50. Fluid-structure interaction solution module; 60. Parameter inversion and correction module; and 70. Early warning visualization module. Detailed Implementation

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

[0049] See Figure 1 The present invention provides a coupled numerical simulation system for formation rainfall-seepage-deformation. This system is configured to integrate multi-source geological monitoring data, multi-field coupled numerical calculation logic, and a physical information-driven parameter inversion algorithm for real-time assessment and prediction of the stability of the Feixianguan Formation in the Wumeng Mountains. The system is based on instructions stored in a computer-readable storage medium and runs on a processor. It connects to a field physical sensor network and a computing terminal via a data bus or communication interface, forming a closed-loop workflow for data acquisition, processing, analysis, and feedback.

[0050] This system is specifically designed for the unique geological structure of the Feixianguan Formation in the Wumeng Mountains. This formation is primarily composed of interbedded soft and hard rocks, with the soft rocks mainly consisting of mudstone or shale, and the hard rocks mainly consisting of limestone or sandstone, all rich in calcareous cement. Given the physicochemical properties of this formation under rainfall and groundwater infiltration, including water-induced disintegration, accelerated chemical dissolution, and interlayer shear slippage, this system does not employ a single mechanical equilibrium calculation. Instead, it constructs a coupled computational architecture that considers shear expansion-induced chemical damage and damage gradient-driven permeability evolution. The system aims to address the challenge of existing geological hazard prediction technologies accurately simulating sudden interlayer slippage and seepage channel connection caused by microstructural changes in this type of formation.

[0051] like Figure 1 The numerical simulation system for coupled formation rainfall-seepage-deformation mainly includes a data acquisition and monitoring module 10, a geological modeling module 20, a constitutive coupling calculation module 30, a seepage channel reconstruction module 40, a fluid-structure interaction solution module 50, a parameter inversion and correction module 60, and an early warning visualization module 70. The modules are connected via data interfaces to achieve one-way data transmission or two-way interaction. The data acquisition and monitoring module 10 is configured to acquire physical quantities of the geological environment; the geological modeling module 20 is configured to provide the geometric model and mesh data required for calculation; the constitutive coupling calculation module 30 and the seepage channel reconstruction module 40 are configured to calculate the dynamic evolution of the material properties of the soil and rock mass; the fluid-structure interaction solution module 50 is configured to calculate the flow field distribution and provide feedback on pore water pressure; the parameter inversion and correction module 60 is configured to correct model parameters based on observation errors; and the early warning visualization module 70 is configured to output the final analysis results.

[0052] Through the coordinated operation of the aforementioned modules, the system can simulate the entire process of how the shear expansion effect caused by interlayer displacement within the Feixianguan Formation accelerates chemical dissolution through a negative pressure mechanism during rainfall infiltration, thereby altering the permeability of the rock mass. The system uses a parameter inversion correction module 60 to compare the measured displacements acquired by the data acquisition and monitoring module 10 with the calculated values, and employs probabilistic statistical methods to correct key physicochemical parameters in the constitutive coupling calculation module 30. This reduces the uncertainty of model predictions and enables early identification and quantitative warning of geological hazards such as landslides.

[0053] See Figure 1 The present invention provides a numerical simulation system for coupled formation rainfall, seepage and deformation, which includes a data acquisition and monitoring module 10, a geological modeling module 20, a constitutive coupling calculation module 30, a seepage channel reconstruction module 40, a fluid-structure coupling solution module 50, a parameter inversion and correction module 60, and an early warning visualization module 70.

[0054] The data acquisition and monitoring module 10 is configured as the system's basic data input terminal. This module connects to a sensor network and geological database distributed in the Feixianguan Formation field in the Wumeng Mountain area via a communication interface. Specifically, the data acquisition and monitoring module 10 is responsible for collecting two types of data: the first type is static geological exploration data, including lithological distribution of rock strata, thickness of interlayers of soft and hard rock, stratum occurrence elements, and laboratory-determined rock physical and mechanical reference parameters obtained from borehole core samples; the second type is dynamic real-time monitoring data, including surface displacement time series obtained through the Global Navigation Satellite System (GNSS), horizontal displacement at different underground depths obtained through a deep inclinometer, groundwater level changes obtained through a pore water pressure gauge, and rainfall intensity and cumulative rainfall data obtained through a rain gauge. The data acquisition and monitoring module 10 performs analog-to-digital conversion and preprocessing on the above analog signals, and transmits the standardized data to the geological modeling module 20 and the parameter inversion and correction module 60.

[0055] The geological modeling module 20 is configured to construct a three-dimensional numerical calculation model reflecting the geological structural characteristics of the target area based on the static geological survey data input from the data acquisition and monitoring module 10. The geological modeling module 20 spatially discretizes the computational domain, generating a finite element mesh or a set of discrete element particles. Considering the unique soft-hard interlayered structure of the Feixianguan Formation, the geological modeling module 20 specifically sets interface elements with independent constitutive properties at the interface between the mudstone soft layer and the limestone or sandstone hard layer. These interface elements are used to capture interlayer shear slip, opening, and sliding deformation behaviors. The three-dimensional geological mesh model with interface elements generated by the geological modeling module 20 is transmitted to the constitutive coupling calculation module 30 and the fluid-structure coupling solution module 50 as the computational domain carrier.

[0056] In this application, the soil-rock mass element adopts an elastoplastic damage constitutive model based on the Mohr-Coulomb criterion. Its stress-strain relationship is expressed as:

[0057] ;

[0058] in, For the effective stress tensor, For the isotropic elastic stiffness tensor, For the total strain tensor, For plastic strain tensor; Represents the total damage variable of the rock mass element at the current moment; symbol " "Represents the double dot product operation of a tensor (the same applies below).

[0059] For the interface elements at the junction of the soft mudstone layer and the hard limestone layer, a thickness-free Cohesive contact model (or Goodman joint element model) is adopted. The stress-displacement relationship is defined as follows:

[0060] ;

[0061] in, and These are the tangential stress and normal stress at the interface, respectively. and These are the relative shear displacement and the normal opening displacement of the interface, respectively. and These are the tangential stiffness and normal stiffness of the interface, respectively. and These are the cohesion and internal friction angle at the interface, respectively. When the shear stress satisfies... At that time, the interface slips, and the stiffness... It degenerates into residual stiffness.

[0062] Constitutive coupling calculation module 30 is configured to calculate the damage evolution state of soil and rock elements under complex stress fields and hydrochemical environments. This module embeds a shear-dissolution strongly coupled algorithm to quantify the interaction between stress damage and chemical damage. Within each calculation time step, constitutive coupling calculation module 30 calculates the plastic shear strain and the associated shear expansion rate based on the current stress state. Through preset chemical kinetic logic, this module identifies the local negative pressure effect caused by shear expansion and calculates the acceleration factor of the chemical dissolution rate accordingly, thereby updating the total damage variables of the soil and rock elements. Constitutive coupling calculation module 30 outputs the updated damage scalar field data to the permeability channel reconstruction module 40.

[0063] The permeability channel reconstruction module 40 is configured to dynamically update the permeability properties of the medium based on the spatial distribution characteristics of rock mass damage. This module receives damage field data from the constitutive coupling calculation module 30 and calculates the gradient vector of the damage scalar field in three-dimensional space. The permeability channel reconstruction module 40 constructs a local coordinate rotation matrix based on the direction of the damage gradient, rotating the principal axis of the permeability tensor to be consistent with the direction of the damage gradient or a preset fracture propagation direction.

[0064] Let the unit normal vector calculated from the damage gradient be... .

[0065] To construct a local coordinate system, we first select an auxiliary vector. (like Minimum, then Otherwise, take unit vectors along other axes, and construct tangential basis vectors through cross product. and .

[0066] Rotation matrix for transforming from global coordinate system to local damage principal axis coordinate system The construction is as follows:

[0067] ;

[0068] This matrix is ​​used to rotate the diagonalized local permeation tensor to the global coordinate system.

[0069] Simultaneously, this module nonlinearly increases the principal value of the permeation tensor based on the magnitude of the damage variable, thereby simulating a directional dominant flow channel in the numerical model. The updated permeation tensor field data is then transmitted to the fluid-structure interaction solution module 50.

[0070] The fluid-structure interaction solution module 50 is configured to solve the fluid dynamics equations for porous media to obtain the pore water pressure distribution. This module calculates the local Reynolds number of the fluid within each grid cell and determines the flow regime based on the comparison between the Reynolds number and a critical threshold. When the Reynolds number is below the threshold:

[0071] This threshold is typically set as the critical Reynolds number. In this embodiment, take The values ​​between (specific values) This threshold is not simply set arbitrarily, but is based on the Forchheimer number. The nonlinear significance is determined by the local Reynolds number. The calculation formula is:

[0072] ;

[0073] in, For fluid density, For seepage velocity, For dynamic viscosity, This refers to the average pore size or fracture aperture of the medium. When the calculated... When the inertial force is negligible compared to the viscous force, it is determined to be linear laminar flow; otherwise, it is determined to be nonlinear turbulent flow.

[0074] This module uses the Darcy's law algorithm for solving the problem; when the Reynolds number is higher than or equal to the threshold, the module automatically switches to the Forchheimer nonlinear seepage algorithm to simulate the high-speed non-Darcy flow effect in the dominant flow channel.

[0075] when In this case, Darcy's Law is applied:

[0076] ;

[0077] when At that time, the Forchheimer nonlinear seepage equation is used:

[0078] ;

[0079] in, For pressure gradient, The permeability of the medium. Non-Darcy inertia coefficient (unit: m) −1 Its empirical formula is usually taken as , here The shape factor, Porosity.

[0080] The pore water pressure field calculated by the fluid-structure interaction solution module 50 is fed back to the constitutive coupling calculation module 30 to update the effective stress, and is also transmitted to the parameter inversion correction module 60.

[0081] The parameter inversion correction module 60 is configured to correct key parameters of the numerical model using field measured data. This module receives measured surface and deep displacement values ​​provided by the data acquisition and monitoring module 10, as well as predicted displacement values ​​output by the fluid-structure interaction solution module 50. The parameter inversion correction module 60 constructs a likelihood function based on physical information constraints, employing a Bayesian inference algorithm or a Markov chain Monte Carlo sampling method.

[0082] Bayesian inference is a parameter estimation method based on probability and statistics. Its core formula is the posterior probability. .in, It is the prior distribution of the parameters (the range of values ​​set based on experience). It is the likelihood function (reflecting the degree of agreement between the model output and the measured data).

[0083] Since the posterior probability distribution is usually extremely complex and cannot be solved analytically, this invention employs the Markov Chain Monte Carlo (MCMC) method for numerical sampling. MCMC constructs a Markov chain whose stationary distribution converges to the target posterior distribution, thereby extracting a large number of samples in the parameter space and using the statistical characteristics of these samples (such as mean and variance) to approximate the optimal solution of the parameters.

[0084] The chemical dissolution rate constant and the mechanical coupling coefficient in the constitutive coupling calculation module 30 are estimated posteriorly by probability. This module outputs the corrected optimal parameter combination to the constitutive coupling calculation module 30 for prediction calculation in the next time step.

[0085] The early warning visualization module 70 is configured to receive and process the system's final calculation results. This module maps the displacement field, stress field, damage field, and pore water pressure field within the strata to a three-dimensional visualization interface. The early warning visualization module 70 has pre-set stability criteria. When the cumulative damage level of a continuous surface exceeds a critical threshold, or when a sudden change in the displacement rate of a key monitoring point occurs, the module automatically generates a landslide disaster early warning signal, identifies potential landslide paths and the estimated failure time, and outputs this information to the user via a display terminal.

[0086] See appendix Figure 2 The present invention provides a numerical simulation method for coupled formation rainfall, seepage and deformation, the method including step S1 initialization and modeling.

[0087] In step S1, the geological modeling module 20 first establishes a three-dimensional numerical calculation model of the Feixianguan Formation in the Wumeng Mountain area. The geological modeling module 20 reads borehole exploration data, topographic data, and stratum attitude data output by the data acquisition and monitoring module 10 through a data interface. Based on the above data, the geological modeling module 20 reconstructs the three-dimensional geometry of the strata and precisely cuts out the spatial distribution areas of soft and hard rock layers through Boolean operations. Considering the structural characteristics of the Feixianguan Formation's interbedded soft and hard rock layers, the geological modeling module 20 sets interface units with independent mechanical properties on the contact surface between the mudstone soft layer and the limestone or sandstone hard layer. These interface units are configured to allow shear displacement and normal opening between layers to simulate various contact behaviors along the bedding planes.

[0088] The interface elements employ contact logic constructed using the penalty method. In the local coordinate system, the constitutive equations of the interface elements are in incremental form:

[0089] ;

[0090] in, This represents the normal relative displacement increment (opening or closing). This represents the tangential relative slip increment.

[0091] To simulate interlaminar damage, tangential stiffness Will follow plastic shear work It evolved through accumulation: ,in For initial stiffness, This is the stiffness degradation parameter.

[0092] The geological modeling module 20 discretizes the constructed geometric model, generating a hybrid mesh containing rock mass units and interlayer interface units. Subsequently, the geological modeling module 20 initializes the material parameters for each computational unit in the model. These parameters include the rock's physical, mechanical, and hydraulic parameters. Physical parameters include dry density, saturated density, and porosity. Mechanical parameters include bulk modulus, shear modulus, internal cohesion, and internal friction angle. Hydraulic parameters include initial matrix permeability, specific storage coefficient, and Biot coefficient. Specifically, for interface units, tangential stiffness, normal stiffness, and interfacial friction coefficient are also assigned. Simultaneously, the geological modeling module 20 sets the initial values ​​of the total damage variable, mechanical damage variable, and chemical damage variable within all computational units to zero, or to the initial background damage values ​​calculated based on in-situ acoustic wave testing.

[0093] After assigning parameter values, the geological modeling module 20 sets the initial boundary conditions of the model. Regarding mechanical boundary conditions, the bottom surface of the model is set as a fully constrained boundary, the sides of the model are set as normal constrained boundaries, and the top surface of the model is set as a free-deformation surface. Regarding hydraulic boundary conditions, the geological modeling module 20 sets constant head boundaries or impermeable boundaries on the sides of the model based on the initial groundwater level monitoring data provided by the data acquisition and monitoring module 10.

[0094] Subsequently, the constitutive coupling calculation module 30 and the fluid-structure interaction solution module 50 perform initial field equilibrium calculations based on the aforementioned mesh and boundary conditions. First, the geostress field is initialized. The system performs elastic calculations under gravity load only until the maximum unbalanced force within the model drops below the preset convergence tolerance, thereby obtaining the initial geostress distribution. The calculated displacement generated at this point is then cleared to eliminate spurious deformations caused by gravity loading. Next, the initial seepage field is calculated. The system solves the steady-state seepage equation to obtain the pore water pressure distribution field corresponding to the initial groundwater level. This initial geostress field and initial pore water pressure field will serve as the starting state for subsequent time steps of rainfall-seepage coupling calculations.

[0095] See appendix Figure 2 The present invention provides a numerical simulation method for coupled formation rainfall, seepage and deformation, which includes step S2: data input and boundary update.

[0096] In step S2, the system enters the time-step cyclic calculation process. First, the data acquisition and monitoring module 10 acquires the external environmental load data and the formation response monitoring data within the current time step. The data acquisition and monitoring module 10 reads the real-time rainfall intensity data recorded by the rain gauge installed on site and transmits the data to the fluid-structure interaction solution module 50. The fluid-structure interaction solution module 50 identifies the upper surface boundary nodes of the three-dimensional numerical model constructed by the geological modeling module 20 and converts the received rainfall intensity values ​​into flow flux boundary conditions. The fluid-structure interaction solution module 50 determines whether the current rainfall intensity exceeds the infiltration capacity of the surface soil. If it does not exceed the infiltration capacity, the rainfall intensity is directly applied as a second type of flow flux boundary condition; if it exceeds the infiltration capacity, the boundary condition is modified to a first type of constant head boundary condition to simulate surface water infiltration, thereby completing the dynamic update of the model's hydraulic boundary.

[0097] Simultaneously, the data acquisition and monitoring module 10 reads the current GNSS displacement reading, the horizontal displacement reading from the deep inclinometer, and the pressure reading from the underground pore water pressure gauge. The data acquisition and monitoring module 10 performs time synchronization processing on the aforementioned multi-source heterogeneous data to ensure that all monitoring data correspond to the same moment in the current simulation time step. The processed monitoring data is transmitted to the parameter inversion and correction module 60. The parameter inversion and correction module 60 constructs an observation vector from the received measured values ​​of surface displacement, deep displacement, and pore water pressure, and stores this observation vector in the system cache as the reference true value for evaluating model prediction bias and performing Bayesian parameter inversion in subsequent steps. Through step S2, the system ensures that the numerical simulation process is always driven by real external meteorological conditions and always has the latest measured data for correcting the model state.

[0098] See appendix Figure 2 The present invention provides a numerical simulation method for coupled formation rainfall, seepage and deformation, the method including step S3 shear-dissolution coupled calculation step.

[0099] In step S3, the constitutive coupling calculation module 30 performs stress update and damage evolution calculations for the soil and rock elements. The constitutive coupling calculation module 30 first receives the strain increment, pore water pressure, and current total damage variable transmitted from the previous time step or initial conditions. This module employs an elastoplastic constitutive algorithm, first calculating the elastic trial stress based on Hooke's law.

[0100] Elastic test stress The calculation is based on the incremental form of the generalized Hooke's law:

[0101] ;

[0102] in, The stress tensor of the previous time step, This represents the total strain increment in the current step. The elastic stiffness matrix is ​​derived from Young's modulus. Compared to Poisson Sure.

[0103] Subsequently, the constitutive coupling calculation module 30 substitutes the elastic test stress into the yield function for discrimination.

[0104] This embodiment uses the Mohr-Coulomb yield criterion as the plasticity criterion, and its yield function... The expression is:

[0105] ;

[0106] in, For average normal stress, This refers to the generalized shear stress (Mises equivalent stress). Lode angle, It is the internal friction angle. For cohesion These are the hardening / softening parameters.

[0107] If the test stress is within the yield surface, the element is determined to be in an elastic state and no new plastic damage occurs; if the test stress exceeds the yield surface, the element is determined to enter a plastic yield state.

[0108] For rock mass elements determined to be yielded, the constitutive coupling calculation module 30 uses a backward Euler return mapping algorithm to pull the stress back to the yield surface.

[0109] when Plasticity correction is required at this time. The backward Euler-backward mapping algorithm is used to solve the nonlinear equation system:

[0110] ;

[0111] in, For the plastic multiplier to be determined, The plastic potential function (usually taken as the non-associated flow rule, with the same form as the yield function but using the dilatation angle) is the plastic potential function. Replace the friction angle The system of equations was solved using Newton's iterative method to obtain the actual stress. .

[0112] The module calculates the increments of plastic shear strain and plastic volumetric strain within that time step. During this process, the constitutive coupling calculation module 30 quantifies the volumetric expansion component caused by shear slip, i.e., shear expansion volumetric strain, based on the non-associated flow law. This module extracts the plastic volumetric strain rate as a key physical quantity characterizing the shear expansion strength for subsequent accelerated calculations of chemical damage.

[0113] The constitutive coupling calculation module 30 then performs chemical damage evolution calculations. This module incorporates a chemical damage acceleration evolution equation based on shear dilatation, which quantitatively describes the enhancing effect of the local negative pressure suction effect generated by the opening of microfractures and the increase in reaction surface area on the chemical dissolution rate during shear dilatation of the rock mass. The constitutive coupling calculation module 30 calculates the chemical damage evolution rate at the current time step using the following formula:

[0114] ;

[0115] in, This represents the derivative of the chemical damage variable with respect to time, i.e., the chemical damage evolution rate. This represents the baseline chemical dissolution rate constant, which is related to the rock mineral composition and the pH value of groundwater. This represents the total damage variable of the rock mass element at the current moment; This represents the mechanocoupler coefficient, used to modulate the acceleration sensitivity of shear expansion to chemical reactions; This represents the plastic volumetric strain rate induced by shear action, i.e., the shear expansion rate.

[0116] The constitutive coupling calculation module 30 calculates the chemical damage increment according to the above formula and linearly superimposes it with the mechanical damage increment calculated by plastic shear strain, thereby updating the total damage variable of the rock mass unit. The updated total damage variable is stored at each integration point and transmitted as an output parameter to the permeability channel reconstruction module 40 to drive subsequent permeability anisotropy evolution calculations. Through step S3, the system realizes the unidirectional strong coupling drive of the mechanical deformation field to the chemical dissolution field, accurately simulating the accelerated dissolution phenomenon caused by shear displacement in the stress concentration zone of the Feixianguan Formation.

[0117] See appendix Figure 2 The present invention provides a numerical simulation method for coupled formation rainfall, seepage and deformation, the method including step S4 dynamic reconstruction of permeability tensor.

[0118] In step S4, the permeability channel reconstruction module 40 performs anisotropic update calculations of rock mass permeability properties. The permeability channel reconstruction module 40 first reads the total damage variable distribution data for the current time step output by the constitutive coupling calculation module 30. The permeability channel reconstruction module 40 employs a numerical differentiation algorithm.

[0119] For a three-dimensional eight-node hexahedral element, the damage gradient at the element center The calculation is performed using the shape function derivative interpolation method:

[0120] ;

[0121] in, For the first The shape function of each node. Local natural coordinates, Global coordinates For the first The damage value at each node (obtained by extrapolation from the integration point).

[0122] The gradient vector of the total damage variable at the center point of each grid cell is calculated in three-dimensional space. This gradient vector indicates the direction of the most drastic change in the damage field, which physically corresponds to the normal direction of the localized shear zone or fracture surface. The permeation channel reconstruction module 40 normalizes this gradient vector to determine the normal basis vector of the local damage coordinate system.

[0123] The permeability channel reconstruction module 40, based on the normal basis vector, constructs a tangential plane basis vector perpendicular to the normal using the Schmidt orthogonalization method, thereby establishing a complete local orthogonal coordinate system. Within this local coordinate system, the permeability channel reconstruction module 40 calculates the permeability coefficient along the tangential direction (dominant flow direction) and along the normal direction of the fracture surface, according to preset damage and permeability evolution rules. The system sets the tangential permeability coefficient to increase nonlinearly exponentially with increasing damage variables to simulate the dominant water-conducting channel effect formed after fracture penetration; while the normal permeability coefficient maintains a lower growth rate or remains at the matrix level.

[0124] Subsequently, the permeation channel reconstruction module 40 constructs a coordinate rotation matrix to transform the diagonalized permeation tensor defined in the local coordinate system back to the global coordinate system, thereby obtaining the effective permeation coefficient matrix in full tensor form. The permeation channel reconstruction module 40 performs permeation tensor rotation and update based on the damage gradient using the following formula:

[0125] ;

[0126] in, This represents the updated global effective penetration tensor; This represents the rotation matrix from the global coordinate system to the local damage coordinate system, which is composed of the direction cosines of the damage gradient vector; Indicates the initial matrix permeability of the rock material; Represents the total damage variable of the current unit; and These represent the permeability growth coefficients parallel to and perpendicular to the damage surface, respectively, used to control the intensity of permeability anisotropy; It represents the power exponent of permeability as a function of damage evolution, and is used to characterize the nonlinearity of the formation of water inrush channels; This represents matrix multiplication. This represents the matrix transpose operation.

[0127] In step S4, the system reconstructs a directional permeability tensor field at each integration point. This updated permeability tensor field reflects the directional seepage channels spontaneously formed within the rock mass due to interlayer slippage and chemical dissolution, and transmits this data to the fluid-structure interaction solution module 50, enabling subsequent flow field calculations to accurately capture high-speed non-Darcy flow phenomena along fault zones or weak interlayers.

[0128] See appendix Figure 2 The present invention provides a numerical simulation method for coupled formation rainfall, seepage and deformation, the method including step S5 nonlinear flow field solution step.

[0129] In step S5, the fluid-structure interaction (FSI) solution module 50 performs iterative solving of the groundwater seepage field based on the updated effective permeability tensor field transmitted by the seepage channel reconstruction module 40. The FSI solution module 50 first determines the fluid flow state within each discrete grid cell. This module reads the fluid seepage velocity vector calculated in the previous iteration step or time step and, combined with the cell characteristic length and fluid viscosity, calculates the local Reynolds number for each cell. The FSI solution module 50 has a built-in critical Reynolds number threshold, which is used to define the transition boundary between laminar and nonlinear turbulent flow within the rock mass. The FSI solution module 50 compares the calculated local Reynolds number with the critical Reynolds number threshold cell by cell.

[0130] When the local Reynolds number of a given element is less than the critical Reynolds number threshold, the fluid-structure interaction solution module 50 determines that the fluid flow within that element conforms to the linear Darcy's law, meaning that the flow velocity and pressure gradient are linearly related. When the local Reynolds number of a given element is greater than or equal to the critical Reynolds number threshold:

[0131] The critical Reynolds number threshold here These are empirical constants determined based on geotechnical hydraulic tests. In practice, the system defaults to these settings. This value can also be obtained by inversion correction using the pressure-flow nonlinear curve from the field pumping test through the parameter inversion correction module 60. Its physical meaning is the critical point when the inertial force term accounts for more than 5%-10% in the momentum equation.

[0132] This indicates that the unit is located within an open fissure or a channel with a high dissolution advantage, where the flow velocity is high and the inertial effect of the fluid flow cannot be ignored. At this point, the fluid-structure interaction solution module 50 automatically switches to the nonlinear Forchheimer seepage algorithm, introducing an inertial drag term to correct the pressure gradient calculation, thereby limiting the permeability at high flow velocities.

[0133] The fluid-structure interaction (FSI) solution module 50 solves the adaptive FSI governing equations that include non-Darcy flow terms. These equations comprehensively describe the evolution of pore water pressure over time and the volumetric coupling effect between the fluid and the solid skeleton. The specific formulas are as follows:

[0134] ;

[0135] in, Indicates the specific storage coefficient of soil and rock media; Indicates pore water pressure; Indicates time; Indicates the effective stress coefficient of Biot; This represents the volumetric strain of the rock mass skeleton, reflecting the cumulative effect of solid deformation on flow field pressure. Represents the divergence operator; This represents the effective permeation tensor input by the permeation channel reconstruction module 40; Indicates the dynamic viscosity of a fluid; This represents the Forchheimer non-Darcy inertia coefficient. It takes a zero value when the Reynolds number is below a threshold and a non-zero positive value when it is above a threshold. Indicates fluid density; The magnitude of the seepage velocity vector; Represents the gravitational acceleration vector; Indicates source and sink terms.

[0136] The fluid-structure interaction solution module 50 uses the Newton-Raphson iterative method to solve the above highly nonlinear control equations.

[0137] The residual scheme after discretization of the fluid-structure interaction control equations is as follows: The Newton-Raphson iterative formula is:

[0138] ;

[0139] in, For the first The pore water pressure vector of the next iteration This is the tangent stiffness matrix (Jacobi matrix). Due to the existence of non-Darcy flow terms, Since the flow rate varies, a nonlinear iterative solution must be used.

[0140] The process continues until the residual norm of the pressure field converges to the preset accuracy. The pore water pressure field data obtained at the current time step is fed back to the constitutive coupling calculation module 30, which is used to update the stress state of the rock mass according to the effective stress principle, thereby completing the strong coupling calculation of the fluid field and the solid deformation field.

[0141] See appendix Figure 2 The present invention provides a numerical simulation method for coupled formation rainfall, seepage and deformation, which includes step S6, Bayesian parameter inversion.

[0142] In step S6, the parameter inversion correction module 60 performs probabilistic inversion and correction calculations of model parameters based on measured data, aiming to reduce the uncertainty of physicochemical parameters that are difficult to directly measure in the deep rock mass. The parameter inversion correction module 60 first determines the target parameter vector to be inverted, which mainly includes the baseline chemical dissolution rate constant and the mechanochemical coupling coefficient defined in the constitutive coupling calculation module 30. The parameter inversion correction module 60 obtains the observation vector for the current time step from the data acquisition and monitoring module 10. This observation vector includes measured values ​​of surface GNSS displacement, deep tilt displacement, and pore water pressure. Simultaneously, this module extracts the calculated values ​​of the model at the corresponding monitoring points to form a simulation vector.

[0143] The parameter inversion correction module 60 constructs a likelihood function to quantify the degree of matching between the simulated vector and the observed vector. This module assumes that the observation error and model structure error follow a multidimensional independent Gaussian distribution and evaluates the confidence level of the current parameter combination by calculating the weighted Euclidean norm of the residual vector. The parameter inversion correction module 60 calculates the likelihood function value for a given parameter combination using the following formula:

[0144] ;

[0145] in, This indicates that, given the model parameter vector Obtain observation data under the conditions The likelihood probability; Indicates the total number of observation data; The covariance matrix represents the observation noise and is used to adjust the impact of different types of monitoring data on the inversion results by weighting. This indicates the corresponding parameter combination. The forward model calculation operator is the simulated output vector obtained by executing steps S3 to S5; This represents the vector transpose operation; This represents the matrix inversion operation.

[0146] Based on the constructed likelihood function and the preset prior distribution of parameters, the parameter inversion correction module 60 uses the Markov chain Monte Carlo (MCMC) sampling algorithm to explore the posterior probability distribution of the parameters. This module employs the Metropolis-Hastings sampling strategy.

[0147] The core of the Metropolis-Hastings algorithm lies in the acceptance rate. The calculation. Let the current parameter state be... From the suggested distribution Extracting candidate parameters The probability of accepting this candidate parameter is:

[0148] ;

[0149] The system generates a random number between [0, 1]. ,like Then let (Accept); otherwise, (Rejected). This mechanism guarantees that the sampling chain eventually converges to the true posterior distribution.

[0150] Candidate samples are generated in the parameter space, and the decision to accept the jump is made based on the ratio of the posterior probability of the candidate sample to that of the current sample. The parameter inversion correction module 60 performs parallel sampling of multiple Markov chains and removes samples in the aging period after reaching a preset number of iterations.

[0151] The parameter inversion correction module 60 performs statistical analysis on the retained sample chain, calculates the expected value or maximum a posteriori estimate of the posterior distribution, and uses it as the corrected optimal physicochemical parameters. The parameter inversion correction module 60 feeds back the corrected baseline chemical dissolution rate constant and the mechanochemical coupling coefficient to the constitutive coupling calculation module 30. Upon receiving the updated parameters, the constitutive coupling calculation module 30 uses the new parameters to perform shear dissolution coupling calculations for the next time step, thus achieving dynamic adaptive calibration of the numerical model as the geological evolution process progresses.

[0152] See appendix Figure 2 The present invention provides a numerical simulation method for coupled formation rainfall, seepage and deformation, the method including step S7 stability judgment and early warning step.

[0153] In step S7, the early warning visualization module 70 performs a stability assessment and disaster risk identification based on the current geological condition. The early warning visualization module 70 first extracts the nodal displacement field, nodal velocity field, and damage variable field at the element integration point from the constitutive coupling calculation module 30 for the current time step. Based on the structural characteristics of the Feixianguan Formation's alternating soft and hard layers, the early warning visualization module 70 automatically identifies potential interlayer slip surface paths and calculates the proportion of the element length with damage variables exceeding the critical failure value along this path to the total length of the slip surface, i.e., the degree of penetration damage. Simultaneously, this module tracks the nodal displacement rate at key monitoring points, which are typically located at the slope toe or geological boundary.

[0154] The early warning visualization module 70 constructs a comprehensive stability judgment criterion, normalizes the calculated penetration damage degree and displacement rate, and then performs weighted fusion to generate the landslide initiation risk index at the current moment. The early warning visualization module 70 calculates the landslide initiation risk index using the following formula:

[0155] ;

[0156] in, This indicates the risk index of landslide initiation; and Let represent the weighting coefficients for penetration damage degree and displacement rate, respectively, and satisfy . ; This indicates the length of the section on the potential sliding surface that has undergone plastic yielding or damage. This represents the total arc length of the potential sliding surface; This represents the instantaneous displacement rate magnitude of key monitoring points; This indicates a preset critical displacement rate threshold, which is determined based on geological history data or soil creep tests.

[0157] The early warning visualization module 70 will calculate the landslide initiation risk index. Compare with the system's preset tiered warning thresholds. When When the value is less than the first threshold, the system determines that the formation is in a stable state; when When the value is between the first and second thresholds, the system outputs a yellow warning signal, indicating the presence of accelerated local deformation; when... When the second threshold is exceeded, the system outputs a red alert and indicates the overall instability area and the expected slip volume. The early warning visualization module 70 renders the aforementioned formation deformation cloud map, plastic zone distribution map, flow velocity vector map, and early warning signal to the user interface.

[0158] After completing the current stability assessment and output, the system checks whether the current simulation time has reached the preset total simulation duration. If not, the system advances the simulation time step by an increment. The program then uses the updated stress field, seepage field, and damage field as the initial state for the next time step, jumps back to step S2 (data input and boundary update), and continues the next round of cyclic calculations. If the total simulation time has been reached, the system terminates the calculation process and automatically generates a comprehensive analysis report.

[0159] Preferably, the present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the described numerical simulation method for coupled formation rainfall, seepage, and deformation. Specific Implementation

[0160] In this embodiment, a highway bedding slope (chainage K120+350) in the Feixianguan Formation of the Wumeng Mountain area was selected as the monitoring and simulation object. The strata of this slope are composed of alternating layers of mudstone (soft rock) and limestone (hard rock), with a dip angle of less than 20 degrees to the slope aspect, indicating a significant risk of interlayer shear displacement.

[0161] See appendix Figure 3 The figure visually illustrates how this invention captures deformation acceleration phenomena that traditional methods cannot predict through coupled computation.

[0162] Initialization and parameter setting (corresponding to step S1): The system first constructs a slope model in the geological modeling module 20 and sets basic parameters based on laboratory tests: initial matrix permeability of mudstone. m / s, baseline chemical dissolution rate constant Force coupling coefficient Penetration growth coefficient power exponent .

[0163] Initial rainfall and linear deformation (corresponding) Figure 3 In the 0-20 hour interval: during the initial stage of rainfall ( Up to 20 hours later, the rainwater had not yet fully infiltrated into the deep, weak interlayers. For example... Figure 3 As shown, the field measured data (square markers) largely coincide with the predicted values ​​of this system (solid line) and the traditional model predicted values ​​(gray dashed line), exhibiting a slow linear increase. At this point, the rock mass is in the elastic deformation stage, and the shear dilatation rate... The chemical damage evolution rate is low, and the damage variable (the dotted line corresponding to the right coordinate axis) remains near the initial background value of 0.05.

[0164] Shear dissolution coupling triggered acceleration (corresponding) Figure 3 (Middle 20-35 hour interval): When the simulation reaches Within hours, rainwater infiltration caused a decrease in the effective stress between layers, and the constitutive coupling calculation module 30 detected plastic shear yielding in the weak zone of the mudstone.

[0165] At this point, the plastic shear expansion coefficient calculated by the system is... The number of cases suddenly increased to 8.0 × 10⁻⁶. −3 The system calculates the chemical damage acceleration factor based on the formula in step S3:

[0166] ;

[0167] This means that due to the negative pressure effect generated by shear expansion, the chemical dissolution rate at this moment is increased to 3.32 times the baseline value.

[0168] As a result, Figure 3 The damage variable curve (dotted line) showed a sharp rise after 24 hours, rapidly climbing from 0.05 to 0.45.

[0169] Permeation channel reconstruction and displacement mutation (corresponding to step S4 and) Figure 3 35 hours later): with damage variables Upon reaching 0.45, the permeability of the permeability channel reconstruction module 40 is updated. The effective permeability increase factor is calculated based on the formula:

[0170] ;

[0171] The interlayer permeability increased by four orders of magnitude instantly, simulating the process of fractures connecting to form water inrush channels.

[0172] This physical process is directly reflected in Figure 3 In the displacement curve: the predicted value (solid line) of this system shows a clear inflection point after 30 hours and rises exponentially, finally reaching a cumulative displacement of about 90 mm after 48 hours, which is highly consistent with the field measured data.

[0173] In contrast, the traditional model's predicted value (gray dashed line) is only 35 mm because it does not consider the coupling of shear dissolution and the dynamic evolution of permeability and is calculated based on the original parameters, thus underestimating the landslide risk.

[0174] Experimental verification and effect comparison:

[0175] To verify the effectiveness of this system, a quantitative comparison was made between the predicted data of this system, the predicted data of the traditional uncoupled model, and the measured data of the field GNSS monitoring station (point G03).

[0176] Traditional model bias analysis: such as Figure 3 As shown, traditional models exhibit significant prediction lag and underestimation of magnitude in the later stages of rainfall. The relative error can reach as high as 60% within hours. This is because traditional models treat soil and rock as materials with constant parameters, and cannot simulate the unique chain reaction of water contact, shearing, and dissolution in the Feixianguan Formation.

[0177] Accuracy Analysis of this System: This system improves accuracy by introducing a parameter inversion correction module 60. Using initial monitoring data at the hourly level Fine-tuning was performed to ensure that the later predicted curve closely followed the measured points. Within hours, the relative error between the predicted value (90mm) and the measured value (92mm) of this system was only 2.2%, successfully providing an early warning of the landslide event.

[0178] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A coupled numerical simulation method for formation rainfall-seepage-deformation, characterized in that, Includes the following steps: A three-dimensional geological model with interbedded soft and hard rocks was established, and the rainfall data collected in the field was converted into hydraulic boundary conditions for the three-dimensional geological model in real time. The soil and rock elements of the three-dimensional geological model adopted an elastoplastic damage constitutive model based on the Mohr-Coulomb criterion. For the interface elements at the junction of the mudstone soft layer and the limestone hard layer, a thickness-free Cohesive contact model was adopted, and its stress-displacement relationship was defined as follows: in, and These are the tangential stress and normal stress at the interface, respectively. and These are the relative shear displacement and the normal opening displacement of the interface, respectively. and These are the tangential stiffness and normal stiffness of the interface, respectively. and These represent the cohesion and internal friction angle at the interface, respectively, when the shear stress satisfies... At that time, the interface slips, and the stiffness... Degenerates into residual stiffness; Calculate the shear dilatation rate of rock mass units in the three-dimensional geological model, thereby determining the chemical damage acceleration factor, and then update the total damage variable of the rock mass units through the chemical damage acceleration factor. The permeation tensor characterizing the dominant flow channel is reconstructed based on the distribution gradient of the total damage variable in three-dimensional space: the unit normal vector calculated from the damage gradient is... Choose an auxiliary vector ,like Minimum, then Otherwise, take unit vectors along other axes and construct tangential basis vectors through cross product. and This allows for the construction of a local coordinate system; a rotation matrix is ​​then used. : The diagonalized local permeability tensor is rotated to the global coordinate system. Then, the fluid-structure interaction equation is solved based on the permeability tensor to obtain the pore water pressure distribution. The local Reynolds number of the fluid is calculated in each grid cell of the three-dimensional geological model. When the local Reynolds number is lower than a preset threshold, the Darcy's law algorithm is called for solving. When the local Reynolds number is higher than or equal to the threshold, the Forchheimer nonlinear seepage algorithm is automatically switched to solve, thereby updating the effective stress and obtaining the calculated displacement field. A likelihood function is constructed based on the deviation between the calculated displacement field and the displacement monitoring data collected on site. The model parameters of the three-dimensional geological model are then corrected by inversion using the likelihood function.

2. The numerical simulation method for coupled formation rainfall-seepage-deformation according to claim 1, characterized in that, Also includes: The modified three-dimensional geological model outputs the rainfall-seepage-deformation coupling results, thereby assessing formation stability and outputting early warning signals.

3. The coupled numerical simulation method for formation rainfall-seepage-deformation according to claim 2, characterized in that, The specific methods for establishing the three-dimensional geological model include: Interpolation algorithms are used to process borehole exploration data to construct stratigraphic interfaces. Geometric cutting operations are then performed on the stratigraphic interfaces to generate soft mudstone layers and hard rock layers. The soft mudstone layers and hard rock layers are then discretized to generate rock mass units in a three-dimensional geological model.

4. The numerical simulation method for coupled formation rainfall-seepage-deformation according to claim 3, characterized in that, The specific method for establishing the three-dimensional geological model further includes: identifying the contact boundary between the soft mudstone layer and the hard rock layer, generating the interface unit at the contact boundary and setting the mechanical parameters of the interface unit, and simulating the shear displacement and slip deformation behavior of the strata along the bedding plane through the interface unit, thereby establishing the three-dimensional geological model.

5. The numerical simulation method for coupled formation rainfall-seepage-deformation according to claim 2, characterized in that, The hydraulic boundary conditions for converting on-site rainfall data into a three-dimensional geological model in real time include: Identify the upper surface boundary nodes of the 3D geological model, extract real-time rainfall intensity from field rainfall data, and compare it with the infiltration capacity of the surface soil at the upper surface boundary nodes; when the real-time rainfall intensity does not exceed the infiltration capacity, apply the real-time rainfall intensity directly as the second type of flow boundary condition; When the real-time rainfall intensity exceeds the infiltration capacity, the hydraulic boundary conditions are modified to the first type of constant head boundary conditions, and the infiltration of surface water is simulated through the first type of constant head boundary conditions.

6. The numerical simulation method for coupled formation rainfall-seepage-deformation according to claim 2, characterized in that, Calculating the shear expansion rate includes: An elastoplastic constitutive algorithm is used to determine whether a rock mass element has entered a plastic yield state. For rock mass elements that have entered a plastic yield state, the volume expansion component caused by shear slip is quantified according to the non-associated flow rule. The plastic volumetric strain rate corresponding to the volume expansion component is extracted and defined as the shear dilatation strain rate.

7. The numerical simulation method for coupled formation rainfall-seepage-deformation according to claim 2, characterized in that, The total damage variable for the rock mass element is updated as follows: A calculation logic for the chemical damage acceleration factor is constructed, which is defined as an exponential function of the natural constant, wherein the exponent of the exponential function is the product of the mechanochemical coupling coefficient and the shear expansion rate.

8. The numerical simulation method for coupled formation rainfall, seepage, and deformation according to claim 7, characterized in that, Calculate the chemical damage evolution rate, which is the product of the baseline chemical dissolution rate constant, a term reflecting the remaining integrity of the rock mass, and the chemical damage acceleration factor; The chemical damage increment is calculated based on the chemical damage evolution rate, and the chemical damage increment is superimposed with the mechanical damage increment to update the total damage variable.

9. The numerical simulation method for coupled formation rainfall, seepage, and deformation according to claim 2, characterized in that, The reconstructed permeation tensor characterizing the dominant flow channel also includes: In the constructed local coordinate system, the tangential permeability coefficient is set to increase non-linearly exponentially with the total damage variable.

10. The numerical simulation method for coupled formation rainfall-seepage-deformation according to claim 2, characterized in that, The process of constructing and correcting the model parameters of a 3D geological model using likelihood functions includes: The measured values ​​of surface displacement and deep displacement are obtained to form an observation vector, and the numerical values ​​of the displacement field at the corresponding locations are extracted to form a simulation vector. Assuming that the observation error follows a multidimensional independent Gaussian distribution, a likelihood function is constructed based on the residual between the observation vector and the simulation vector to quantify the degree of matching.

11. The numerical simulation method for coupled formation rainfall, seepage, and deformation according to claim 10, characterized in that, The construction and inversion correction of the model parameters of the three-dimensional geological model through the likelihood function also includes: generating candidate samples based on the likelihood function using the Markov chain Monte Carlo sampling algorithm, and determining the corrected model parameters through statistical analysis. The model parameters include the baseline chemical dissolution rate constant and the mechanochemical coupling coefficient.

12. The numerical simulation method for coupled formation rainfall, seepage, and deformation according to claim 2, characterized in that, The assessment of formation stability and output of early warning signals includes: Identify potential interlayer sliding surface paths and calculate the proportion of the length of the element whose total damage variable on the interlayer sliding surface path exceeds the critical failure value to the total length of the sliding surface as the penetration damage degree. The displacement rates of key monitoring points are obtained and normalized. The penetration damage degree and the normalized displacement rate are weighted and fused to generate a landslide initiation risk index. When the landslide initiation risk index exceeds a preset threshold, a warning signal containing the potential landslide path and the expected time of failure is generated.

13. A coupled numerical simulation system for formation rainfall, seepage, and deformation, characterized in that, The method for performing the coupled numerical simulation of formation rainfall-seepage-deformation as described in any one of claims 1 to 12 includes: The data acquisition and monitoring module is used to obtain physical quantity inputs of the geological environment, and to collect static geological survey data and dynamic real-time monitoring data. The geological modeling module is used to construct a three-dimensional geological grid model containing interface units based on the data input from the data acquisition and monitoring module. The constitutive coupling calculation module is used to calculate the shear dilatancy rate and chemical damage acceleration factor of rock mass elements, and update the total damage variable. The permeation channel reconstruction module is used to reconstruct the permeation tensor characterizing the dominant flow channel based on the distribution gradient of the total damage variable in three-dimensional space. The fluid-structure interaction solution module is used to solve the fluid-structure interaction equation based on the permeation tensor characterizing the dominant flow channel to obtain the pore water pressure distribution; The parameter inversion correction module is used to construct a likelihood function using field displacement monitoring data and perform inversion correction on the parameters in the constitutive coupling calculation module. The early warning visualization module is used to receive calculation results and generate landslide disaster early warning signals based on stability criteria.

14. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the program is executed by the processor, it implements the numerical simulation method for coupled formation rainfall, seepage and deformation as described in any one of claims 1-12.

Images