Soil-stone medium digital reconstruction and twinborn simulation method

Abstract

The invention relates to the technical field of soil-rock admixtures, and provides a soil-rock medium digital reconstruction and twinborn simulation method, gravel-containing unsaturated soil comprises gravel and unsaturated soil, the gravel part is modeled by adopting a discontinuous deformation analysis (DDA) method, and the movement, contact and interaction of the gravel part are simulated; the unsaturated soil part is discretized by adopting a smoothed particle hydrodynamics (SPH) method, a mixed medium theory is combined, an unsaturated soil moisture motion equation is introduced, moisture migration and pore water pressure change in a soil body are described, and hydraulic response characteristics of the gravel-containing unsaturated soil are reflected by reasonably setting hydraulic boundary conditions; the interaction force between the gravel and the particles is transmitted through the contact module, and the contact module dynamically marks the particles near the block as impermeable particles in the calculation process so as to prevent water phase leakage. According to the device and the method, the seepage-deformation coupling effect simulation of the gravel-containing unsaturated soil can be better realized.

Description

Technical Field

[0001] This invention relates to the field of soil-rock mixture technology, and more specifically, to a method for digital reconstruction and twin simulation of soil-rock media. Background Technology

[0002] Gravelly unsaturated soil, a typical heterogeneous geological body in mountainous areas, is regulated by complex hydrogeological processes in its formation and evolution. From a genetic perspective, gravelly unsaturated soil is formed under the influence of seepage, driven by rainfall infiltration and groundwater level fluctuations, resulting in a multiphase system consisting of a solid phase (gravel and soil), a liquid phase (pore water), and a gaseous phase (pore gas). Under unsteady seepage, the physical and mechanical properties of gravelly unsaturated soil, such as matrix suction and shear strength, deteriorate significantly. Specifically, increased pore water pressure leads to a decrease in effective stress, directly weakening the shear strength of slopes and ultimately inducing landslides, posing a serious threat to life, property, and infrastructure. Therefore, revealing the water-mechanical coupling failure mechanism of gravelly unsaturated soil under seepage is crucial for slope disaster early warning and prevention.

[0003] Numerical simulation is an important tool for studying seepage-induced landslide geological hazards. However, the simulation of seepage in gravelly unsaturated soil faces multiple challenges: significant block effects lead to the coexistence of mechanical properties of discrete and continuous media; complex pore structures result in highly nonlinear seepage paths; strong coupling between water migration and solid skeleton deformation in unsaturated states, as well as the dynamic evolution of large deformation behavior after instability, all pose severe challenges to traditional numerical methods. While continuous medium methods, such as the finite element method (FEM), can handle fluid-structure interaction problems, their assumption of an equivalent continuous medium makes it difficult to accurately describe the strong heterogeneity and discreteness between gravel and soil, as well as the complex contact, collision, and large rotation behaviors among gravel. Their applicability is particularly limited under conditions of high gravel content or large deformation failure. While the Discrete Element Method (DEM) excels at simulating the motion and contact of discrete particles, accurately simulating unsaturated seepage processes involving massive fine-grained soil requires generating a huge number of particle elements and coupling them with computational fluid dynamics (CFD), resulting in extremely high computational costs. This makes it difficult to meet the needs of analysis at the practical engineering scale, and modeling key characteristics such as matrix suction in unsaturated seepage is also quite complex.

[0004] The coupling of Discontinuous Deformation Analysis (DDA) and Smooth Particle Hydrodynamics (SPH) (DDA-SPH) provides an efficient solution for the interaction between irregularly shaped blocks and discrete particles. DDA, based on contact and kinematic constraints between discrete blocks (such as rock), implicitly solves the displacement and deformation of the block system, making it suitable for large deformation and instability problems in discontinuous media such as jointed rock masses. SPH, a meshless Lagrange particle method, directly solves the hydrodynamic equations by approximating particle properties and their spatial derivatives through kernel functions, excelling at simulating rheological problems with large deformations and free surfaces, such as water and soil movement. Currently, DDA-SPH has been successfully applied to two-phase coupled problems such as soil-structure interaction, landslides in soil-rock mixtures, and tunnel water inrush. However, this method has significant limitations in numerically simulating the seepage-deformation coupling process in gravelly unsaturated soils. Such problems involve multiphase interactions between rock blocks (solid phase), soil particles (solid phase), liquid water (liquid phase), and pore gas (gas phase). Existing frameworks struggle to handle matrix suction effects in unsaturated regions, gas-liquid capillary interactions, and strong nonlinear coupling mechanisms between solid deformation. Summary of the Invention

[0005] The present invention provides a method for digital reconstruction and twin simulation of soil-rock media, which can solve the problem of rainfall-induced soil-rock mixed landslides.

[0006] According to a method for digital reconstruction and twin simulation of soil-rock media based on the present invention, the gravel-bearing unsaturated soil includes gravel and unsaturated soil. The gravel part is modeled using the discontinuous deformation analysis method (DDA) to simulate its motion, contact and interaction. The unsaturated soil part is discretized using the smoothed particle hydrodynamics method (SPH). The method combines the mixed media theory and introduces the water movement equation of the unsaturated soil to describe the water migration and pore water pressure changes inside the soil. The hydraulic response characteristics of the gravel-bearing unsaturated soil are reflected by reasonably setting the hydraulic boundary conditions. The interaction force between gravel and particles is transmitted through a contact module. During the calculation process, the contact module dynamically marks particles near the block as non-permeable particles to prevent water phase leakage.

[0007] As a preferred option, the unsaturated soil portion is constructed using a single-layer SPH method to create an unsaturated uniform medium, thereby achieving water-soil mixing. The single-layer SPH method endows each SPH particle with the ability to characterize the macroscopic state of a locally unsaturated soil unit by introducing porosity. n With saturation S r Two state variables are implemented, so a single particle can simultaneously represent a mixed medium consisting of soil particles, pore water, and pore gas. Let the unit volume of unsaturated soil be V According to the mixed-media theory, there is ,in V s , V w and V a Corresponding to soil phases s Liquid phase w Harmony a Volume; density of unsaturated soil The expression is obtained by weighting the density of each component: In the formula, r α The inherent density of each component; n α Volume percentage of each component V α / V And satisfy ; α Indicates the medium type; ignoring the gas phase density, expanding and refining the above equation, we obtain the overall density of the mixture as: in, n Porosity; S r Saturation; r s For soil phase density, r w The density is given by the aqueous phase; based on Terzaghi's effective stress principle, the total stress tensor of unsaturated soil is further derived. s The format is as follows: in, s α Let be the stress characteristic quantities of each phase; neglecting the influence of pore gas pressure, expanding and refining the above equation, we get: in, n Porosity; S r Saturation; p w Pore ​​water pressure; d ij This is the Kronecker tensor, used to describe the isotropic components of the stress tensor.

[0008] As a preferred method, when simulating the dynamic process of rainfall infiltration or groundwater level changes, the migration of water in unsaturated soil strictly follows the governing equations; solving the governing equations drives the water flux between particles, which in turn triggers the dynamic evolution of local particle saturation and pore water pressure; in the calculation, the particles need to be rebalanced based on the continuity equation and momentum equation to strictly maintain the overall mass conservation and momentum conservation of the system.

[0009] As a preferred option, the mass conservation principle for the rebalancing calculation is as follows: Taking a unit volume of unsaturated soil, assuming the material derivatives of the solid, fluid, and air phases are zero, the mass balance equations for each phase are derived. Defining the material derivative of the soil skeleton, the general form of the mass balance equations for each phase of unsaturated soil is as follows:

[0010] in, α Indicates the media type; v α for α The absolute velocity of the phase medium; ∂ (•) / ∂t This is a local time term, describing the change in material mass over time; The convection term reflects the mass transport effect caused by phase migration; D s (•) / Dt This represents the rate of change of a physical quantity with respect to the total fluid particles. D s Indicates other relative solid phases s Find the mass derivative; neglecting the mass of the gas phase and considering only the mass conservation of the soil and water phases, the mass balance equations for the soil and water phases are established as follows: Assuming the soil phase is incompressible, i.e., its density... r s The mass balance equation of the soil phase, which does not change with time or space, is transformed into the form of a porosity change rate: The absolute velocity of the water phase is decomposed into the sum of the skeleton velocity and the relative velocity; the volumetric water content of the unsaturated soil is known. i = nS r The water phase mass balance equation is transformed into:

[0011] In the formula, v ws =v w - v s This is the velocity of the water phase relative to the soil skeleton, also known as the Darcy velocity; further assuming the water phase has a uniform density, i.e. Using the material derivative form with the soil skeleton as a reference, the equilibrium equation is transformed into: Under quasi-static conditions, the change in pore water pressure and density satisfy a linear relationship, combined with the bulk modulus of water. K w The association is established as follows: The relationship between saturation and pore water pressure is as follows: In the formula, p c For the matrix suction of unsaturated soil, p c = p a - p w By combining the soil and water phase mass balance equations, the water density-pressure relationship, and the saturation-pressure relationship, and eliminating intermediate variables, the pore water pressure governing equation is derived as follows: In the formula, ∂θ / ∂p c It is a characteristic parameter of the soil-water characteristic curve of unsaturated soil, reflecting the volumetric water content. i Sensitivity to matrix suction.

[0012] As a preferred approach, the momentum conservation in the rebalancing calculation is as follows: As a multiphase medium, unsaturated soil experiences water and gas migration, with water and gas flowing within the pores of the soil particle skeleton. This process is accompanied by complex interaction forces. The momentum conservation equation can accurately describe these interactions. In the mixed medium, each phase must obey the linear momentum conservation equation to ensure the transfer and balance of momentum during deformation and flow. In the mixed medium, neglecting the gas phase effect, the momentum conservation equations for each phase are as follows: in, b It is a volume force; for α Harmony β The viscous drag force generated by the contact between phase media, calculated according to Darcy's law, is expressed as follows: In the formula, k It is the permeability coefficient; g It is the acceleration due to gravity; based on the above equation, the linear momentum balance equation for the water phase is derived: Assuming quasi-static conditions, we have dv w / dt =0; further decompose the stress tensor of the aqueous phase into s w ≈- p w I, where I is the identity matrix; thus, the Darcy velocity of pore water motion... Represented as: Known pw = ghh ,in h For pressure head; and ∇ z =-1, z Let the position be the water head, with upward being positive; thus: Therefore, the Darcy speed above further transforms into: Ignoring the effect of pore gas pressure, i.e. yes =0, then the matrix suction pc =- pw =- ghh , h For pressure head, yw Let the natural unit weight of water be denoted; neglecting soil skeleton deformation and water phase compressibility, the pore water pressure governing equation degenerates into the Richards equation, as shown below: in, C ( h ( ) is the specific water capacity. C ( h )=( ∂θ / ∂h This reflects the sensitivity of the unsaturated soil volumetric water content to changes in pressure head; furthermore, the overall mixture system must also satisfy overall momentum balance, as shown in the following equation: The Richards equation is essentially a combination of Darcy's law and the law of conservation of mass, used to describe the governing equations of water movement in unsaturated porous media; its core function is to describe the flow process of water in unsaturated soil pores driven by both matrix suction and gravity.

[0013] As a preferred approach, unsaturated seepage in soil manifests as a transient flow process in anisotropic porous media. When solving using the Smooth Particle Hydrodynamics (SPH) method, the governing equations need to be transformed into a system of ordinary differential equations. By selecting an appropriate kernel function, the characteristic parameters of adjacent particles are weighted and averaged to discretize, thus achieving a numerical characterization of the dynamic behavior of the mixed media. According to the basic principles of SPH, without considering energy exchange between materials, the particle density needs to be dynamically updated. r ,speed v ,stress s and positional state parameters; for unsaturated porous media, porosity n Pressure head h and pore water pressure pw It is a key hydraulic parameter affecting its physical and mechanical properties; by using the particle approximation of SPH, the rate of change of porosity is transformed into a discrete form in Cartesian coordinates, as shown in the following equation: particle j Represented as particles i Supports contributing particles within the domain; mj Represents particles j quality pj Represents particles j The density; via , yes Represents particles i , j exist α The velocity component in the direction; We Represents the kernel function; xia Represents particles i exist α The coordinate components of the direction; When pore air pressure is neglected, pore water pressure and pressure head satisfy the following conditions: pw = ρwgh Therefore, only the pressure head needs to be adjusted. h The evolution equation is solved discretely; the pressure head change rate in the above equation is expanded into a discrete form of SPH; the harmonic mean is used to evaluate the permeability coefficient. k Symmetry processing is used to improve computational stability, resulting in... dh / dt The discrete form can be written as follows:

[0014] In the formula, or To prevent the denominator from being zero, the term is a minus sign. ; Yes As a source-sink term, when considering water inflow and infiltration, the source is used to replenish the water content; H Given the total head; the second-order precision leapfrog method (LF) is used to update the state variables. The step-by-step update format of the LF algorithm is as follows:

[0015]

[0016] n and h These are soil particle porosity and hydraulic pressure head, respectively. t The number of steps at that time, t +1 / 2 is the previous step. t To the next step t The median of +1.

[0017] As a preferred option, the equations for the movement of water in unsaturated soil include soil phase constitutive equations and hydraulic control equations. The soil constitutive equation is as follows: Since gravel in gravelly unsaturated soils is non-absorbent, the unsaturated characteristic of the mixture is mainly manifested in the matrix suction of the unsaturated soil. According to Bishop's effective stress principle for unsaturated soils, when pore air pressure is not considered, the effective stress tensor of unsaturated soils expressed by the SPH method... s ´ is represented as: Where χ is the Bishop effective stress parameter; the effective stress tensor will serve as the fundamental variable describing the stress-strain relationship of the soil skeleton; the Drucker-Prager yield criterion-related flow law is used to describe the elastoplastic deformation of unsaturated soil; considering the large deformation problem of the soil, the Jaumann stress rate is introduced into the constitutive equation, and the stress-strain relationship is expressed as: in, K Bulk modulus; G Shear modulus; J 2 is the second invariant of the deviatoric stress tensor; It is the deviatoric strain rate tensor; This is the shear stress tensor; ds / dt The time derivative of stress; The time derivative of the strain; The time derivative of the rotation; superscript α ,β and c A tensor dummy index; and It is a calculation coefficient; f The internal friction angle of the soil; particles i strain rate tensor and spin rate tensor They are represented as follows: Among them, particles i The velocity gradient is expressed as:

[0018] The hydraulic governing equations are as follows: Water transport and storage in gravelly unsaturated soils occur only in the unsaturated soil phase. The VG model is used to characterize the soil-water properties of unsaturated soils, and its volumetric water content... i With pressure head h The relationship is as follows: In the formula, l , g These are the fitting parameters for soil-related characteristic curves; i s This represents the saturated volumetric water content. i r The residual volumetric water content is represented by the hydraulic conductivity curve, which is a predictive model derived from the VG model combined with Mualem pore distribution theory, used to describe the hydraulic conductivity coefficient in unsaturated soil. k With effective saturation S r e The relationship between the changes is expressed as follows: In the formula, k sat is the saturated permeability coefficient of the soil.

[0019] As a preferred approach, considering the characteristics of the research object and working conditions, a seepage model for gravelly unsaturated soil based on the DDA-SPH method is constructed, specifically defining the following types of hydraulic boundary conditions: (1) Non-drainage boundary conditions Undrained boundary conditions refer to applying a constraint that the normal fluid flux is zero at the model boundary. This means that water is not allowed to flow perpendicularly through the interface into or out of the computational domain, but flow along the boundary tangentially is permitted. Based on the application scenario and location, undrained boundary conditions are divided into the following two categories: 1) Category I non-drainage boundary This type of boundary is used to simulate truncated boundaries to maintain stable pore water pressure at the boundary; to achieve the condition of zero normal flux, three layers of fixed virtual particles are set outside the model boundary, and their pressure head is set accordingly. hj and position of water head zj Setting and integration domain center particle i The corresponding values ​​are the same, that is hj = hi , zj = zi This ensures that the total head of the boundary virtual particles is equal to that of the internal real particles, i.e. Hj = Hello Under this setting, the key head difference term controlling the normal flux in the SPH discrete scheme is always zero, i.e. Hj - Hello =0, thus strictly satisfying the condition. qn The physical requirement of =0; 2) Second type of non-drainage boundary This type of boundary is used to define the contact interface between the soil and the impermeable rock block inside, to prevent water from seeping into the rock block region. A boundary kernel function truncation method is used to limit the hydraulic interaction between particles and the rock block. This treatment may lead to incomplete particle support domains near the boundary, causing numerical leakage and mass non-conservation. Therefore, it is necessary to further refine the boundary conditions for particles near the rock block. i The kernel function was renormalized, and the corrected kernel function is as follows:

[0020] We Represents the original kernel function. N For particles i Contributing particles within the integration domain j The number; by using the normalized modified kernel function Wijcor This compensates for the loss of kernel function weights caused by boundary effects, suppresses numerical instability and non-conservation of mass, and ensures that there is no water leakage at the boundaries of internal rock blocks. (2) Free drainage boundary conditions The bottom boundary of the gravelly unsaturated soil model is defined as a free drainage boundary to simulate the natural drainage process; this boundary follows Darcy's law, and the virtual particle pressure head outside the boundary is constant at atmospheric pressure. hj =0, and satisfies the unidirectional drainage condition: when the total internal head Hello The total head of virtual particles above the boundary Hj hour, Hello > Hj Water is freely expelled; when Hello < Hj At that time, normal flux qn=0; To achieve this constraint, three layers of fixed virtual particles are set outside the boundary; the pressure head of the virtual particles is set to 0. hj =0, position head zj Then it is precisely set according to the actual terrain elevation; in the SPH discrete format, the internal real particles i With boundary virtual particles j Total head difference Δ Hi =( hi + zi )- zj Automatically adjust flux direction: when Δ Hi When Δ > 0, the SPH diffusion term generates a positive flux to achieve drainage; when Δ Hi When ≤0, the gradient of the kernel function is ∇ We The dot product with the position vector makes the flux term approach zero, thus strictly satisfying the condition. qn One-way drainage condition ≥0; (3) Boundary conditions for surface rainfall infiltration The rainfall infiltration boundary in a gravelly unsaturated soil model was implemented using the direct source term method to simulate the effect of rainfall on soil moisture content. During the simulation, rainfall was treated as a source term and directly incorporated into the discrete form of the Richards equations. For surface particles, the rainfall source term... Yes Represented as: In the formula, qrain Rainfall intensity; Hey For particles i Surface area; Vi This represents the volume of soil particles; to prevent rainfall intensity from exceeding the soil's infiltration capacity, an actual infiltration constraint is introduced, defining the actual infiltration volume. qactual Rainfall intensity qrain The smaller of the soil permeability coefficient and the total soil permeability coefficient; the specific expression is:

[0021] in who The unsaturated permeability coefficient, ksat The saturated permeability coefficient, hi For pressure head; (4) Constant head boundary To accurately characterize the hydraulic impact of the initial groundwater level on the system, a fixed constant head boundary is set at the model's side boundary: based on the spatial distribution of the initial groundwater level, a hydraulic reference surface is constructed by assigning a constant head value to the boundary virtual particles; simultaneously, a two-way drainage and seepage boundary condition is established, using the hydraulic gradient to drive the two-way seepage flow exchange between the external water body and the model interior, and finally, the system seepage balance is achieved by relying on the dynamic feedback mechanism of the boundary head difference and seepage flow, thus fully simulating the steady-state control effect of the groundwater level; when the dynamic response of the groundwater level needs to be considered, the constant head is replaced with a periodic time-varying head boundary to realize the continuous evolution simulation of the time-varying effects of groundwater.

[0022] The beneficial effects of this invention are as follows: This invention first uses the single-layer SPH method to construct an unsaturated uniform medium to achieve water-soil mixing; then, through the coupling technology of SPH and DDA, it simulates the interaction between the water-soil mixed porous medium and the rocks, and finally realizes the simulation of the seepage-deformation coupling effect of gravelly unsaturated soil.

[0023] This invention, based on the DDA-SPH coupled framework, introduces an unsaturated porous media seepage model into the SPH module, incorporating the physical processes characterizing water-air two-phase migration. DDA is used to simulate gravel, with its boundary defined as an impermeable boundary, establishing a multiphase interaction criterion for gravel-soil-water in unsaturated gravel-containing soil. By introducing seepage-deformation control equations, a coupled solution for "seepage force-driven → fine particle migration → gravel movement" is ultimately achieved. Simulations based on unsaturated soil-stone column infiltration tests validate the model's reliability and are further applied to a rainfall-induced mixed slope landslide case, successfully reproducing the entire process from initial instability to kinetic deposition, revealing a chain-like catastrophe mechanism of "seepage front expansion → soil and rock mass strength decay → slope instability and failure." This method provides a new theoretical tool for the mechanism research and risk prevention of soil-rock mixed landslides. Attached Figure Description

[0024] Figure 1 This is a flowchart illustrating a method for digital reconstruction and twin simulation of soil-rock media in an embodiment; Figure 2 This is a schematic diagram of DDA-SPH coupling and particle labeling in the embodiment; Figure 3 This is a schematic diagram of the DDA-SPH gravel-containing unsaturated soil seepage model considering rainfall in the embodiments. Detailed Implementation

[0025] To further understand the content of this invention, a detailed description of the invention will be provided in conjunction with the accompanying drawings and embodiments. It should be understood that the embodiments are merely illustrative and not limiting of the invention.

[0026] Example Gravel-bearing unsaturated soil is mainly composed of four phases: gravel, soil, water, and gas. Clearly, existing single numerical simulation techniques are insufficient to accurately simulate all four phases simultaneously, primarily due to the complex multiphase coupling computational requirements. This embodiment assumes that the gravel is an ideal medium that is impermeable to water and air, in which case water and gas are mainly contained within the soil, constituting unsaturated soil. Based on this, this embodiment proposes a digital reconstruction and twin simulation method for soil-rock media. Gravel-bearing unsaturated soil includes both gravel and unsaturated soil. The gravel portion is modeled using Discontinuous Deformation Analysis (DDA) to simulate its motion, contact, and interaction. The unsaturated soil portion is discretized using Smooth Particle Hydrodynamics (SPH). Combining mixed media theory and introducing the unsaturated soil water movement equation, the method describes the water migration and pore water pressure changes within the soil. By appropriately setting hydraulic boundary conditions, the hydraulic response characteristics of gravel-bearing unsaturated soil are reflected. The interaction forces between the gravel and the particles are transmitted through contact modules, such as... Figure 2 As shown, the contact module dynamically marks particles near the block as non-permeable particles during the calculation process to prevent water phase leakage.

[0027] DDA method Discontinuous variable analysis (DDA) aims to solve the problems of large deformation and instability in discontinuous media such as jointed rock masses. DDA is a block system dynamics method based on implicit time integration and the principle of minimum potential energy. For two-dimensional DDA, the basic principle is to discretize the rock mass into independent polygonal blocks, each block... i centroid displacement vector D i Defined by 6 degrees of freedom: in, u0 and v0 Is the center of mass in x and y The translational displacement in the direction; r0 is the rotational displacement of the block about its centroid; ex , ey It is a positive strain; c xy represents shear strain. The total potential energy Π of the system consists of the block strain energy, contact force potential energy, and work done by external forces. By minimizing the total potential energy using equation (3), the overall equilibrium state is obtained: Where K is the global stiffness matrix; F is the load vector. For a block system containing n blocks, its governing equilibrium equations can be expressed in matrix form as follows:

[0028] Where Kij is the stiffness submatrix; Fi is the generalized load vector on block i. The centroid displacement D of the block is solved by "open-closed iteration". i Applying the complete first-order displacement approximation, any point within block i ( x , y displacement () u , v It can be achieved through the displacement transformation matrix T i calculate: in, x 0、 y 0 represents the coordinates of the block's centroid.

[0029] SPH method The core breakthrough of Smooth Particle Hydrodynamics (SPH) lies in breaking free from mesh constraints, becoming the first meshless Lagrangian particle method, and successfully solving the mesh distortion defects of traditional hydrodynamic methods in problems such as ultra-large deformation and free surface fragmentation.

[0030] The SPH method is based on the principles of the kernel approximation and the particle approximation. Its core idea is to discretize the continuous medium into mass-carrying components. m ,density r ,pressure p The mathematical essence of a swarm of particles with similar physical properties is based on the use of smooth kernel functions. W By applying a weighted average to neighborhood particles, the field function and its derivative at any location in space can be reconstructed. Kernel function W It is a weighted function used to smooth the physical quantities in the region surrounding a particle. (Arbitrary continuous function) f ( r ) and its derivative It can be represented as: in, Oh Represent the problem domain. h Indicates the definition of a kernel function W The smooth length of the support field, ∇ r Indicates relative position r The gradient. Smooth length. h This is an important parameter that determines the influence range of the kernel function, thus affecting the resolution and accuracy of the simulation. The support domain of a particle refers to the region of influence within a certain distance when calculating particle interactions. h Other particles within the support domain. In the standard SPH framework, the function values ​​and derivatives of particles are obtained by interpolating the values ​​of other particles within the support domain, as follows: Among them, subscript i and j Indicates particle index, N Represents particles i The total number of particles within the supported domain. mj and pj They are particles j The mass and density. For soil, the differential terms in the Navier-Stokes equations (NS) can be transformed into particle summation form as shown below, achieving discretization of the dynamic equations.

[0031] Among them, m , r and v These represent mass, density, and velocity, respectively. s Indicates stress; F This indicates that the force is external, such as boundary contact force. g Π represents volume force; Π is the artificial viscosity term, introduced by Monaghan to suppress numerical oscillations.

[0032] The DDA-SPH coupling method is a powerful numerical tool for solving complex problems involving the intense, discontinuous interaction between discrete solid blocks and continuous fluids (or fluid-like substances), which are difficult to handle effectively using traditional single methods. It has unique advantages and promising applications in simulating hydraulic structure failure, geological hazards, and impact penetration. With the improvement of computing power and the continuous refinement of coupling algorithms, its application scope will become increasingly broad.

[0033] Rebalancing of characteristic parameters To accurately characterize the complex interactions between soil particles, pore water, and pore gas in unsaturated soil, this study proposes a single-layer SPH method based on mixture theory.

[0034] The single-layer SPH method endows each SPH particle with the ability to characterize the macroscopic state of a locally unsaturated soil unit by introducing porosity. n With saturation S r Two state variables are implemented, so a single particle can simultaneously represent a mixed medium consisting of soil particles, pore water, and pore gas. Let the unit volume of unsaturated soil be V According to the mixed-media theory, there is ,in V s , Vw and V a Corresponding to soil phases s Liquid phase w Harmony a The volume of unsaturated soil; the density of unsaturated soil is obtained by weighting the densities of its components, expressed as: In the formula, r α The inherent density of each component; n α Volume percentage of each component V α / V And satisfy ; α Indicates the medium type; ignoring the gas phase density, expanding and refining the above equation, we obtain the overall density of the mixture as: in, n Porosity; S r Saturation; r s For soil phase density, r w Let be the density of the aqueous phase; based on Terzaghi's effective stress principle, the total stress tensor of unsaturated soil is further derived, in the following form: in, s α Let be the stress characteristic quantities of each phase; neglecting the influence of pore gas pressure, expanding and refining the above equation, we get: in, n Porosity; S r Saturation; p w Pore ​​water pressure; d ij This is the Kronecker tensor, used to describe the isotropic components of the stress tensor.

[0035] In simulating the dynamic process of rainfall infiltration or groundwater level fluctuation, the migration of water in unsaturated soil strictly follows the governing equations; the solution of the governing equations drives the water flux (mass exchange) between particles, thereby triggering the dynamic evolution of local particle saturation and pore water pressure; in the calculation, the particles need to be rebalanced based on the continuity equation and momentum equation to strictly maintain the overall mass conservation (water migration) and momentum conservation (mechanical response) of the system. In this embodiment, the numerical model has been appropriately simplified, and the derivation of the governing equations follows the following basic assumptions: (1) neglecting gas phase pressure and mass ( p a =0, r a =0); (2) Soil and water phases are incompressible and their inherent density is constant; (3) The fluid is inviscid; (4) There is no seepage inside the rock; (5) Darcy's law of seepage is satisfied.

[0036] The mass conservation principle for rebalancing is as follows: Taking a unit volume of unsaturated soil, assuming the material derivatives of the solid, fluid, and air phases are zero, the mass balance equations for each phase are derived. Defining the material derivative of the soil skeleton, the general form of the mass balance equations for each phase of unsaturated soil is as follows:

[0037] in, α Indicates the media type; v α for α The absolute velocity of the phase medium; ∂ (•) / ∂t This is a local time term, describing the change in material mass over time; v s •∇(•) represents the convection term, reflecting the mass transport effect caused by phase migration; D s (•) / Dt This represents the rate of change of a physical quantity with respect to the total fluid particles. D s Indicates other relative solid phases s Find the mass derivative; neglecting the mass of the gas phase and considering only the mass conservation of the soil and water phases, the mass balance equations for the soil and water phases are established as follows: Assuming the soil phase is incompressible, i.e., its density... r s The mass balance equation of the soil phase, which does not change with time or space, is transformed into the form of a porosity change rate:

[0038] The absolute velocity of the water phase is decomposed into the sum of the skeleton velocity and the relative velocity; the volumetric water content of the unsaturated soil is known. i = nS r The water phase mass balance equation is transformed into: In the formula, v ws = v w - v s The velocity of the water phase relative to the soil skeleton, also known as the Darcy velocity; further assuming that the water phase has a uniform density, i.e.,  r w =0, using the material derivative form with the soil skeleton as a reference, the equilibrium equation is transformed into: Under quasi-static conditions, the change in pore water pressure and density satisfy a linear relationship, combined with the bulk modulus of water. K w The association is established as follows: The relationship between saturation and pore water pressure is as follows: In the formula, p c For the matrix suction of unsaturated soil, p c = p a - p w By combining the soil and water phase mass balance equations, the water density-pressure relationship, and the saturation-pressure relationship, and eliminating intermediate variables, the pore water pressure governing equation is derived as follows: In the formula, ∂θ / ∂p c It is a characteristic parameter of the soil-water characteristic curve of unsaturated soil, reflecting the volumetric water content. i Sensitivity to matrix suction.

[0039] The momentum conservation in the rebalancing calculation is as follows: As a multiphase medium, unsaturated soil experiences water and gas migration, with water and gas flowing within the pores of the soil particle skeleton. This process involves complex interactions. The momentum conservation equation precisely describes these interactions, such as the resistance encountered by water and gas flowing through the pores by the soil particle skeleton and the drag force generated between the two phases. In the mixed medium, each phase must obey the linear momentum conservation equation to ensure the transfer and balance of momentum during deformation and flow. In the mixed medium, neglecting the gas phase effect, each phase (soil particle phase)... s Aqueous phase w The momentum conservation equation of () is in the form of: in, b It is a volume force; for α Harmony β The viscous drag force generated by the contact between phase media, calculated according to Darcy's law, is expressed as follows: In the formula, k It is the permeability coefficient; g It is the acceleration due to gravity; based on the above equation, the linear momentum balance equation for the water phase is derived: Assuming quasi-static conditions, we have dv w / dt =0; further decompose the stress tensor of the aqueous phase into s w ≈- p w I, where I is the identity matrix; therefore, the Darcy velocity of pore water motion is expressed as: Known pw = ghh ,in h For pressure head; and ∇ z =-1, z Let the position be the water head, with upward being positive; then we can obtain: Therefore, the Darcy speed above can be further transformed into: Ignoring the effect of pore gas pressure, i.e. yes =0, then the matrix suction pc =- pw =- ghh , h For pressure head, cw is the natural unit weight of water; neglecting soil skeleton deformation and water phase compressibility (∇vs=0, Kw→∞), the pore water pressure governing equation degenerates into the Richards equation, as shown below: in, C ( h ( ) is the specific water capacity. This reflects the sensitivity of the unsaturated soil volumetric water content to changes in pressure head; furthermore, the overall mixture system must also satisfy overall momentum balance, as shown in the following equation: The Richards equation is essentially a combination of Darcy's law and the law of conservation of mass, used to describe the governing equations for water movement in unsaturated porous media. Its core function is to describe the flow process of water in unsaturated soil pores driven by both matrix suction (or negative pore water pressure) and gravity. It is a key theoretical tool for predicting the spatiotemporal distribution of water in the unsaturated zone under the influence of rainfall infiltration, evaporation, and groundwater level fluctuations.

[0040] Discrete form of SPH for hydraulic characteristic parameters Unsaturated seepage in soil manifests as a transient flow process in anisotropic porous media. When solving this process using the Smooth Particle Hydrodynamics (SPH) method, the governing equations need to be transformed into a system of ordinary differential equations. By selecting a suitable kernel function, the characteristic parameters of adjacent particles are weighted and averaged to discretize, thus achieving a numerical characterization of the dynamic behavior of the mixed media. According to the basic principles of SPH, without considering energy exchange between materials, the particle density needs to be dynamically updated. r ,speed v ,stress s and state parameters such as location; for unsaturated porous media, porosity n Pressure head h and pore water pressure pw It is a key hydraulic parameter affecting its physical and mechanical properties; through the particle approximation of SPH, the rate of change of porosity (dn / dt) is transformed into a discrete form in Cartesian coordinates, as shown in the following equation: particle j Represented as particles i Supports contributing particles within the domain; mj Represents particles j quality pj Represents particles j The density; via , yes Represents particles i , j exist α The velocity component in the direction; We Represents the kernel function; xia Represents particles i exist The coordinate components of the direction; When pore air pressure is neglected, pore water pressure and pressure head satisfy the following conditions: pw = ρwgh Therefore, only the pressure head needs to be adjusted. h The evolution equation is solved discretely; the pressure head change rate in the above equation is expanded into a discrete form of SPH; to improve the computational stability of the discrete equation, a harmonic average is used for the permeability coefficient. k Symmetry processing is used to improve computational stability, resulting in... dh / dt The discrete form can be written as follows:

[0041] In the formula, or To prevent the denominator from being zero, the term is a minus sign. ; Yes As a source-sink term, when considering water inflow and infiltration, the source is used to replenish the water content; H Given the total head; the second-order precision leapfrog method (LF) is used to update the state variables. The step-by-step update format of the LF algorithm is as follows:

[0042]

[0043] n and h These are soil particle porosity and hydraulic pressure head, respectively. t The number of steps at that time, t +1 / 2 is the previous step. t To the next step t The median of +1.

[0044] Hydraulic and mechanical characteristics of gravelly unsaturated soil The equations governing water movement in unsaturated soil include the soil constitutive equation and the hydraulic control equation; the soil constitutive equation is specifically as follows: Since gravel in gravelly unsaturated soils is non-absorbent, the unsaturated characteristic of the mixture is mainly manifested in the matrix suction of the unsaturated soil. According to Bishop's principle of effective stress in unsaturated soils, when pore air pressure is not considered, the effective stress tensor σ´ of unsaturated soils expressed by the SPH method can be expressed as: in, x These are Bishop effective stress parameters, typically x ≈ SrThe effective stress tensor will serve as the fundamental variable describing the stress-strain relationship of the soil skeleton; the Drucker-Prager yield criterion-related flow law will be used to describe the elastoplastic deformation of unsaturated soil; considering the large deformation problem of the soil, the Jaumann stress rate will be introduced into the constitutive equation, and the stress-strain relationship will be expressed as: in, K Bulk modulus; G Shear modulus; J 2 is the second invariant of the deviatoric stress tensor; It is the deviatoric strain rate tensor; This is the shear stress tensor; ds / dt The time derivative of stress; The time derivative of the strain; The time derivative of the rotation; superscript α , β and c A tensor dummy index; and It is a calculation coefficient; f The internal friction angle of the soil; particles i strain rate tensor and spin rate tensor They are represented as follows: Among them, particles i The velocity gradient is expressed as: The hydraulic governing equations are as follows: Water transport and storage in gravelly unsaturated soils occur only in the unsaturated soil phase. The VG model is used to characterize the soil-water characteristic curve of unsaturated soils, and its water content... i With pressure head h The relationship is as follows: In the formula, l , g These are the fitting parameters for soil-related characteristic curves; θs This represents the saturated volumetric water content. θr The residual volumetric water content is represented by the hydraulic conductivity curve, which is a predictive model derived from the VG model combined with Mualem pore distribution theory, used to describe the hydraulic conductivity coefficient in unsaturated soil. k With effective saturation Sre The relationship between the changes is expressed as follows: In the formula, ksat The saturated permeability coefficient of the soil; Sre It is the effective saturation; l This is a parameter related to pore connectivity.

[0045] Based on the characteristics of the research object and working conditions, a seepage model for gravelly unsaturated soil was constructed using the DDA-SPH method, such as... Figure 3 As shown, the following types of hydraulic boundary conditions are specifically defined: (1) Non-drainage boundary conditions Undrained boundary conditions refer to applying a zero normal fluid flux to the model boundary. qn The constraint ≡0) means that water is not allowed to flow vertically through the interface into or out of the computational domain, but flow along the boundary tangentially is allowed. Based on the application scenario and location, this paper classifies undrained boundary conditions into the following two categories: 1) Category I non-drainage boundary This type of boundary is used to simulate truncated boundaries to maintain stable boundary pore water pressure, such as Figure 3 As shown in (a); to achieve the condition of zero normal flux, three layers of fixed virtual particles are set outside the model boundary, and their pressure head is set accordingly. hj and position of water head zj Setting and integration domain center particle i The corresponding values ​​are the same, that is hj = hi , zj = zi This ensures that the total head of the boundary virtual particles is equal to that of the internal real particles. Hj = Hello Under this setting, the key head difference term controlling the normal flux in the SPH discrete scheme ( Hj - Hello The value is always zero, thus strictly satisfying the condition. qn The physical requirement of =0; 2) Second type of non-drainage boundary This type of boundary is used to define the contact interface between the soil and the impermeable rock blocks inside, to prevent water from seeping into the rock block area, such as... Figure 3 As shown in (c); given the complex shape of the rock, to avoid the high computational cost of the mirror particle method, this study adopts the boundary kernel function truncation method to limit the hydraulic interaction between particles and the rock; this treatment may lead to incomplete particle support domains near the boundary, causing numerical leakage and mass non-conservation; therefore, it is necessary to further refine the particle support domains near the rock. iThe kernel function was renormalized, and the corrected kernel function is as follows:

[0046] We Represents the original kernel function. N For particles i Contributing particles within the integration domain j The number; by using the normalized modified kernel function Wijcor It can compensate for the loss of kernel function weights caused by boundary effects, suppress numerical instability and mass non-conservation, and ensure that there is no water leakage at the boundary of internal rock blocks. (2) Free drainage boundary conditions The bottom boundary of the gravelly unsaturated soil model is defined as a free drainage boundary, such as... Figure 3 As shown in section (d), this is used to simulate the natural drainage process; the boundary follows Darcy's law, and the virtual particle pressure head outside the boundary is constant at atmospheric pressure, i.e. hj =0, and satisfies the unidirectional drainage condition: when the total internal head Hello The total head of virtual particles above the boundary Hj hour, Hello > Hj Water is freely expelled; when Hello < Hj At that time, normal flux qn =0; To achieve this physical constraint, three layers of fixed virtual particles are set outside the boundary; the pressure head of the virtual particles is set to 0. hj =0, position head zj Then it is precisely set according to the actual terrain elevation; in the SPH discrete format, the internal real particles i With boundary virtual particles j Total head difference Δ Hi =( hi + zi )- zj Automatically adjust flux direction: when Δ Hi When Δ > 0, the SPH diffusion term generates a positive flux to achieve drainage; when Δ Hi When ≤0, the gradient of the kernel function is  We The dot product with the position vector makes the flux term approach zero, thus strictly satisfying the condition. qn One-way drainage condition ≥0; (3) Surface rainfall infiltration boundary Rainfall infiltration is a flux boundary condition, such as Figure 3 As shown in section (b), the rainfall infiltration boundary of the gravelly unsaturated soil model is realized using the direct source term method to simulate the effect of rainfall on soil moisture content. During the simulation, rainfall is treated as a source term and directly added to the discrete form of the Richards equations. For surface particles, the rainfall source term... Yes Represented as: In the formula, qrain Rainfall intensity; Hey For particles i Surface area ( A =△ d ×1, △ d (Initial particle spacing); Vi Represents the volume of soil particles ( Vi = my / p , my and p These are the overall particle mass and overall density, respectively. To prevent rainfall intensity from exceeding the soil's infiltration capacity, an actual infiltration constraint is introduced, defining the actual infiltration amount. qactual Rainfall intensity qrain The smaller of the soil permeability coefficient and the total soil permeability coefficient; the specific expression is: in who The unsaturated permeability coefficient, ksat The saturated permeability coefficient, hi For pressure head; (4) Constant head boundary To accurately characterize the hydraulic impact of the initial groundwater level on the system, a fixed constant head boundary is set at the model's side boundary: based on the spatial distribution of the initial groundwater level, a hydraulic reference surface is constructed by assigning a constant head value to the boundary virtual particles; simultaneously, a two-way drainage and seepage boundary condition is established, using the hydraulic gradient to drive the two-way seepage flow exchange between the external water body and the model interior, and finally, the system seepage balance is achieved by relying on the dynamic feedback mechanism of the boundary head difference and seepage flow, thus fully simulating the steady-state control effect of the groundwater level; when the dynamic response of the groundwater level needs to be considered, the constant head is replaced with a periodic time-varying head boundary to realize the continuous evolution simulation of the time-varying effects of groundwater.

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

Claims

1. A method for digital reconstruction and twin simulation of soil-rock media, characterized in that: Gravel-bearing unsaturated soil includes gravel and unsaturated soil. The gravel part is modeled using discontinuous deformation analysis (DDA) to simulate its motion, contact and interaction. The unsaturated soil part is discretized using smoothed particle hydrodynamics (SPH). Combined with mixed medium theory and the unsaturated soil water movement equation, the water migration and pore water pressure changes inside the soil are described. By reasonably setting hydraulic boundary conditions, the hydraulic response characteristics of gravel-bearing unsaturated soil are reflected. The interaction force between gravel and particles is transmitted through a contact module. During the calculation process, the contact module dynamically marks particles near the block as non-permeable particles to prevent water phase leakage.

2. The method for digital reconstruction and twin simulation of soil-rock media according to claim 1, characterized in that: The unsaturated soil portion was constructed using a single-layer SPH method to create an unsaturated uniform medium, achieving water-soil mixing. The single-layer SPH method endows each SPH particle with the ability to characterize the macroscopic state of a locally unsaturated soil unit by introducing porosity. n With saturation S r Two state variables are implemented, so a single particle can simultaneously represent a mixed medium consisting of soil particles, pore water, and pore gas. Let the unit volume of unsaturated soil be V According to the mixed-media theory, there is ,in V s , V w and V a Corresponding to soil phases s Liquid phase w Harmony a Volume; density of unsaturated soil The expression is obtained by weighting the density of each component: ； In the formula, ρ α The inherent density of each component; n α Volume percentage of each component V α / V And satisfy ; α Indicates the medium type; ignoring the gas phase density, expanding and refining the above equation, we obtain the overall density of the mixture as: ； in, n Porosity; S r Saturation; ρ s For soil phase density, ρ w The density is given by the aqueous phase; based on Terzaghi's effective stress principle, the total stress tensor of unsaturated soil is further derived. σ The format is as follows: ； in, σ α Let be the stress characteristic quantities of each phase; neglecting the influence of pore gas pressure, expanding and refining the above equation, we get: ； in, n Porosity; S r Saturation; p w Pore ​​water pressure; δ ij This is the Kronecker tensor, used to describe the isotropic components of the stress tensor.

3. The method for digital reconstruction and twin simulation of soil-rock media according to claim 2, characterized in that: When simulating the dynamic process of rainfall infiltration or groundwater level changes, the migration of water in unsaturated soil strictly follows the governing equations. Solving the governing equations drives the water flux between particles, which in turn triggers the dynamic evolution of local particle saturation and pore water pressure. In the calculation, the particles need to be rebalanced based on the continuity equation and momentum equation to strictly maintain the overall mass conservation and momentum conservation of the system.

4. The method for digital reconstruction and twin simulation of soil-rock media according to claim 3, characterized in that: The mass conservation principle for rebalancing is as follows: Taking a unit volume of unsaturated soil, assuming the material derivatives of the solid, fluid, and air phases are zero, the mass balance equations for each phase are derived. Defining the material derivative of the soil skeleton, the general form of the mass balance equations for each phase of unsaturated soil is as follows: ； ； in, α Indicates the media type; v α for α The absolute velocity of the phase medium; ∂ (•) / ∂t This is a local time term, describing the change in material mass over time; The convection term reflects the mass transport effect caused by phase migration; D s (•) / Dt This represents the rate of change of a physical quantity with respect to the total fluid particles. D s Indicates other relative solid phases s Find the mass derivative; neglecting the mass of the gas phase and considering only the mass conservation of the soil and water phases, the mass balance equations for the soil and water phases are established as follows: ； ； Assuming the soil phase is incompressible, i.e., its density... ρ s The mass balance equation of the soil phase, which does not change with time or space, is transformed into the form of a porosity change rate: ； The absolute velocity of the water phase is decomposed into the sum of the skeleton velocity and the relative velocity; the volumetric water content of the unsaturated soil is known. θ = nS r The water phase mass balance equation is transformed into: ； In the formula, v ws = v w - v s The velocity of the water phase relative to the soil skeleton, also known as the Darcy velocity; further assuming that the water phase has a uniform density, i.e.,  ρ w =0, using the material derivative form with the soil skeleton as a reference, the equilibrium equation is transformed into: ； Under quasi-static conditions, the change in pore water pressure and density satisfy a linear relationship, combined with the bulk modulus of water. K w The association is established as follows: ； ； The relationship between saturation and pore water pressure is as follows: ； In the formula, p c For the matrix suction of unsaturated soil, p c = p a - p w By combining the soil and water phase mass balance equations, the water density-pressure relationship, and the saturation-pressure relationship, and eliminating intermediate variables, the pore water pressure governing equation is derived as follows: ； In the formula, ∂θ / ∂p c It is a characteristic parameter of the soil-water characteristic curve of unsaturated soil, reflecting the volumetric water content. θ Sensitivity to matrix suction.

5. The method for digital reconstruction and twin simulation of soil-rock media according to claim 4, characterized in that: The momentum conservation in the rebalancing calculation is as follows: As a multiphase medium, unsaturated soil experiences water and gas migration, with water and gas flowing within the pores of the soil particle skeleton. This process is accompanied by complex interaction forces. The momentum conservation equation can accurately describe these interactions. In the mixed medium, each phase must obey the linear momentum conservation equation to ensure the transfer and balance of momentum during deformation and flow. In the mixed medium, neglecting the gas phase effect, the momentum conservation equations for each phase are as follows: ； in, b It is a volume force; for α Harmony β The viscous drag force generated by the contact between phase media, calculated according to Darcy's law, is expressed as follows: ； In the formula, k It is the permeability coefficient; g It is the acceleration due to gravity; based on the above equation, the linear momentum balance equation for the water phase is derived: ； Assuming quasi-static conditions, we have dv w / dt =0; further decompose the stress tensor of the aqueous phase into σ w ≈- p w I, where I is the identity matrix; thus, the Darcy velocity of pore water motion... Represented as: ； Known pw = γwh ,in h For pressure head; and ∇ z =-1, z Let the position be the water head, with upward being positive; thus: ； Therefore, the Darcy speed above further transforms into: ； Ignoring the effect of pore gas pressure, i.e. pa =0, then the matrix suction PC =- pw =- γwh , h For pressure head, γw Let the natural unit weight of water be denoted; neglecting soil skeleton deformation and water phase compressibility, the pore water pressure governing equation degenerates into the Richards equation, as shown below: ； in, C ( h ( ) is the specific water capacity. C ( h )=( ∂θ / ∂h This reflects the sensitivity of the unsaturated soil volumetric water content to changes in pressure head; furthermore, the overall mixture system must also satisfy overall momentum balance, as shown in the following equation: ； The Richards equation is essentially a combination of Darcy's law and the law of conservation of mass, used to describe the governing equations of water movement in unsaturated porous media; its core function is to describe the flow process of water in unsaturated soil pores driven by both matrix suction and gravity.

6. The method for digital reconstruction and twin simulation of soil-rock media according to claim 5, characterized in that: Unsaturated seepage in soil manifests as a transient flow process in anisotropic porous media. When solving this process using the Smooth Particle Hydrodynamics (SPH) method, the governing equations need to be transformed into a system of ordinary differential equations. By selecting a suitable kernel function, the characteristic parameters of adjacent particles are weighted and averaged to discretize, thus achieving a numerical characterization of the dynamic behavior of the mixed media. According to the basic principles of SPH, without considering energy exchange between materials, the particle density needs to be dynamically updated. ρ ,speed v ,stress σ and positional state parameters; for unsaturated porous media, porosity n Pressure head h and pore water pressure pw It is a key hydraulic parameter affecting its physical and mechanical properties; by using the particle approximation of SPH, the rate of change of porosity is transformed into a discrete form in Cartesian coordinates, as shown in the following equation: ； particle j Represented as particles i Supports contributing particles within the domain; mj Represents particles j quality ρj Represents particles j The density; vi α , vjα Represents particles i , j exist α The velocity component in the direction; Wij Represents the kernel function; xiα Represents particles i exist α The coordinate components of the direction; When pore air pressure is neglected, pore water pressure and pressure head satisfy the following conditions: pw = ρwgh Therefore, only the pressure head needs to be adjusted. h The evolution equations are solved discretely. Expand the pressure head change rate in the above equation into a discrete form of SPH; use the harmonic mean to adjust the permeability coefficient. k Symmetry processing is used to improve computational stability, resulting in... dh / dt The discrete form can be written as follows: ； In the formula, η To prevent the denominator from being zero, the term is a minus sign. ; Si As a source-sink term, when considering water inflow and infiltration, the source is used to replenish the water content; H Given the total head; the second-order precision leapfrog method (LF) is used to update the state variables. The step-by-step update format of the LF algorithm is as follows: ； ； n and h These are soil particle porosity and hydraulic pressure head, respectively. t The number of steps at that time, t +1 / 2 is the previous step. t To the next step t The median of +1.

7. The numerical simulation method for discontinuous deformation of gravelly unsaturated soil considering seepage effects according to claim 6, characterized in that: The equations for water movement in unsaturated soil include soil phase constitutive equations and hydraulic control equations; The soil constitutive equation is as follows: Since gravel in gravelly unsaturated soils is non-absorbent, the unsaturated characteristic of the mixture is mainly manifested in the matrix suction of the unsaturated soil. According to Bishop's effective stress principle for unsaturated soils, when pore air pressure is not considered, the effective stress tensor of unsaturated soils expressed by the SPH method... σ ´ is represented as: ； Where χ is the Bishop effective stress parameter; the effective stress tensor will serve as the fundamental variable describing the stress-strain relationship of the soil skeleton; the Drucker-Prager yield criterion-related flow law is used to describe the elastoplastic deformation of unsaturated soil; considering the large deformation problem of the soil, the Jaumann stress rate is introduced into the constitutive equation, and the stress-strain relationship is expressed as: ； ； ； in, K Bulk modulus; G Shear modulus; J 2 is the second invariant of the deviatoric stress tensor; It is the deviatoric strain rate tensor; This is the shear stress tensor; dσ / dt The time derivative of stress; The time derivative of the strain; The time derivative of the rotation; superscript α , β and γ A tensor dummy index; and It is a calculation coefficient; φ The internal friction angle of the soil; particles i strain rate tensor and spin rate tensor They are represented as follows: ； ； Among them, particles i The velocity gradient is expressed as: ； The hydraulic governing equations are as follows: Water transport and storage in gravelly unsaturated soils occur only in the unsaturated soil phase. The VG model is used to characterize the soil-water properties of unsaturated soils, and its volumetric water content... θ With pressure head h The relationship is as follows: ； In the formula, λ , ζ These are the fitting parameters for soil-related characteristic curves; θ s This represents the saturated volumetric water content. θ r The residual volumetric water content is represented by the hydraulic conductivity curve, which is a predictive model derived from the VG model combined with Mualem pore distribution theory, used to describe the hydraulic conductivity coefficient in unsaturated soil. k With effective saturation S r e The relationship between the changes is expressed as follows: ； ； In the formula, k sat is the saturated permeability coefficient of the soil.

8. The method for digital reconstruction and twin simulation of soil-rock media according to claim 7, characterized in that: Based on the characteristics of the research object and working conditions, a seepage model for gravelly unsaturated soil is constructed using the DDA-SPH method, and the following types of hydraulic boundary conditions are specifically defined: (1) Non-drainage boundary conditions Undrained boundary conditions refer to applying a constraint that the normal fluid flux is zero at the model boundary. This means that water is not allowed to flow perpendicularly through the interface into or out of the computational domain, but flow along the boundary tangentially is permitted. Based on the application scenario and location, undrained boundary conditions are divided into the following two categories: 1) Category I non-drainage boundary This type of boundary is used to simulate truncated boundaries to maintain stable pore water pressure at the boundary; to achieve the condition of zero normal flux, three layers of fixed virtual particles are set outside the model boundary, and their pressure head is set accordingly. hj and position of water head zj Setting and integration domain center particle i The corresponding values ​​are the same, that is hj = hi , zj = zi This ensures that the total head of the boundary virtual particles is equal to that of the internal real particles, i.e. Hj = Hi Under this setting, the key head difference term controlling the normal flux in the SPH discrete scheme is always zero, i.e. Hj - Hi =0, thus strictly satisfying the condition. qn The physical requirement of =0; 2) Second type of non-drainage boundary This type of boundary is used to define the contact interface between the soil and the impermeable rock block inside, to prevent water from seeping into the rock block region. A boundary kernel function truncation method is used to limit the hydraulic interaction between particles and the rock block. This treatment may lead to incomplete particle support domains near the boundary, causing numerical leakage and mass non-conservation. Therefore, it is necessary to further refine the boundary conditions for particles near the rock block. i The kernel function was renormalized, and the corrected kernel function is as follows: ； Wij Represents the original kernel function. N For particles i Contributing particles within the integration domain j The number; by using the normalized modified kernel function Wijcor This compensates for the loss of kernel function weights caused by boundary effects, suppresses numerical instability and non-conservation of mass, and ensures that there is no water leakage at the boundaries of internal rock blocks. (2) Free drainage boundary conditions The bottom boundary of the gravelly unsaturated soil model is defined as a free drainage boundary to simulate the natural drainage process; this boundary follows Darcy's law, and the virtual particle pressure head outside the boundary is constant at atmospheric pressure. hj =0, and satisfies the unidirectional drainage condition: when the total internal head Hi The total head of virtual particles above the boundary Hj hour, Hi > Hj Water is freely expelled; when Hi < Hj At that time, normal flux qn =0; To achieve this constraint, three layers of fixed virtual particles are set outside the boundary; the pressure head of the virtual particles is set to 0. hj =0, position head zj Then it is precisely set according to the actual terrain elevation; in the SPH discrete format, the internal real particles i With boundary virtual particles j Total head difference Δ Hij =( hi + zi )- zj Automatically adjust flux direction: when Δ Hij When Δ > 0, the SPH diffusion term generates a positive flux to achieve drainage; when Δ Hij When ≤0, the gradient of the kernel function is ∇ Wij The dot product with the position vector makes the flux term approach zero, thus strictly satisfying the condition. qn One-way drainage condition ≥0; (3) Boundary conditions for surface rainfall infiltration The rainfall infiltration boundary in a gravelly unsaturated soil model was implemented using the direct source term method to simulate the effect of rainfall on soil moisture content. During the simulation, rainfall was treated as a source term and directly incorporated into the discrete form of the Richards equations. For surface particles, the rainfall source term... Si Represented as: ； In the formula, qrain Rainfall intensity; Ai For particles i Surface area; Vi Represents the volume of soil particles; To prevent rainfall intensity from exceeding the soil's infiltration capacity, an actual infiltration constraint is introduced, defining the actual infiltration rate. qactual Rainfall intensity qrain And the smaller value of the soil permeability coefficient; The specific expression is: ； in ki The unsaturated permeability coefficient, ksat The saturated permeability coefficient, hi For pressure head; (4) Constant head boundary To accurately characterize the hydraulic impact of the initial groundwater level on the system, a fixed constant head boundary is set at the model's side boundary: based on the spatial distribution of the initial groundwater level, a hydraulic reference surface is constructed by assigning a constant head value to the boundary virtual particles; simultaneously, a two-way drainage and seepage boundary condition is established, using the hydraulic gradient to drive the two-way seepage flow exchange between the external water body and the model interior, and finally, the system seepage balance is achieved by relying on the dynamic feedback mechanism of the boundary head difference and seepage flow, thus fully simulating the steady-state control effect of the groundwater level; when the dynamic response of the groundwater level needs to be considered, the constant head is replaced with a periodic time-varying head boundary to realize the continuous evolution simulation of the time-varying effects of groundwater.

Images