Simulation method and system of ground fissure under pre-existing fracture geological conditions coupled with FEM-PD

Abstract

This invention discloses a method and system for simulating ground fissures under pre-existing fracture geological conditions using a coupled FEM-PD model, relating to the field of numerical simulation technology for ground subsidence and ground fissures. The method includes the following steps: constructing a finite element-peripheral dynamics model and initializing its parameters; statically dividing the geological computational domain into a finite element region, a periphery dynamics region, and a coupling transition region; applying pore water pressure changes as an external load to the finite element region; obtaining the global displacement field and stress field by solving the finite element governing equations; applying the global displacement field as a displacement boundary condition unidirectionally to the boundary of the periphery dynamics region; simulating the motion process of material points using an adaptive dynamic relaxation method; and calculating local damage by using the Mohr-Coulomb failure criterion to determine the bond fracture between two material points. This invention innovatively couples the finite element model with the periphery dynamics model, accurately simulating both regional continuity and local discontinuity deformation targets.

Description

Technical Field

[0001] This invention belongs to the field of numerical simulation technology of ground subsidence and ground fissures, and particularly relates to a method and system for simulating ground fissures under pre-existing fault geological conditions coupled with FEM-PD. Background Technology

[0002] Uneven ground subsidence caused by excessive groundwater extraction is often accompanied by ground fissures, seriously threatening engineering safety and regional sustainable development. Effectively simulating the entire process of ground subsidence disaster transforming into ground fissures is of great significance for disaster prevention and control. Under the condition of pre-existing fracture structures in the rock or soil mass (hereinafter referred to as pre-existing fractures), the distribution of ground subsidence will be more uneven, and the development of ground fissures will be significantly controlled by the pre-existing fractures, tending to extend along them or concentrate at their ends. This poses a significant challenge to accurate simulation.

[0003] Numerical simulation is the primary method for studying this problem. Currently, methods such as the Finite Element Method (FEM), based on the continuum assumption, suffer from limitations in handling large deformations and spontaneous fractures, including difficulty in characterizing discontinuous evolution processes, simulation accuracy dependence on mesh generation, and inability to spontaneously handle crack initiation. Peridynamics (PD) describes the mechanical behavior of materials through spatial integral equations, naturally simulating crack initiation and propagation, but its computational cost is extremely high. To balance efficiency and the ability to simulate fracture, existing research has coupled the Finite Element Method with PD, primarily through dynamic transformation coupling and sub-model coupling.

[0004] For dynamic transformation coupling, this method initially uses an all-FEM mesh. During the calculation, FEM elements that meet the failure criterion are dynamically replaced with PD material points. Although this method can adapt to the crack path, its transformation process depends on stress triggering during the calculation. This means that the initial stage of the simulation cannot perform targeted and refined modeling for known geological weak points or pre-existing fracture structures, making it difficult for the simulation to reflect the guiding role of known geological structures on crack evolution.

[0005] The aforementioned dynamic transformation coupling framework suffers from a significant drawback in balancing computational stability and geological determinism. Its core logic relies on dynamically predicting crack initiation based on stress criteria, but it exhibits considerable limitations when dealing with geological models possessing clearly defined pre-existing fractures. First, its "trigger-based" transformation mechanism involves cumbersome finite element mesh topology reconstruction and real-time mapping of near-field dynamic material points, introducing unnecessary computational redundancy. Second, the transient process of element failure transforming into material points easily induces numerical oscillations in the physical field and energy balance inaccuracies, leading to decreased computational stability and the introduction of numerical uncertainty. More critically, because this method lacks pre-existing information such as fracture location and strike, it cannot accurately characterize the guiding role of known weak zones in the evolution of ground fissures, resulting in low simulation accuracy and poor convergence when simulating complex geological hazards controlled by deterministic structures.

[0006] For sub-model coupling, this method embeds the PD region as a "sub-model" into the FEM master model. First, a global FEM calculation is performed, and then the FEM displacement results are applied as boundary conditions to the PD sub-model for local fracture analysis. While this method is conceptually clear, it typically involves complex bidirectional data transfer or iteration to reflect the impact of fracture on overall stiffness, resulting in a cumbersome computational process.

[0007] The aforementioned sub-model coupling method suffers from shortcomings, including inefficiency and a mismatch in physical mechanisms. Its complex bidirectional iterative coupling mechanism (such as using element deletion to feedback stiffness reduction) aims to simulate the impact of a single brittle crack on the structural bearing capacity. However, the ground settlement-ground fissure problem is essentially a coupled evolution of continuous deformation and discontinuous rupture in large-scale soil and rock masses driven by pore water pressure loads. Its core is the transmission and coordination of loads between continuous and discontinuous media, rather than the weakening of structural stiffness by a single crack. Existing sub-model methods require tedious boundary data interpolation and stiffness matrix reorganization iterations at each time step, attempting to capture the nonlinear feedback of local element failures on the overall structure in real time. However, in large-scale ground settlement simulations, the key lies in solving the stable transmission of macroscopic seepage field loads between continuous grids and discontinuous lattices, rather than the local stiffness oscillations caused by microscopic element deletion. Therefore, the frequent cross-domain data mapping of the sub-model method introduces significant algorithmic redundancy, leading to low computational efficiency when dealing with large-scale field problems. Meanwhile, the common bond elongation rate is often used as the failure criterion in the PD region, which fails to optimize the shear failure mechanism of the soil and rock mass, resulting in insufficient simulation accuracy.

[0008] Therefore, there is an urgent need for a dedicated FEM-PD numerical simulation framework that can fully integrate prior geological information such as existing fractures, possess an efficient interface load transfer mechanism, and deeply integrate the shear failure mechanism of rock and soil with the fluid-structure interaction effect, so as to achieve accurate characterization of the initiation and propagation process of ground fissures under the background of large-scale ground subsidence.

[0009] To address the limitations of existing finite element-peripheral dynamics (FEM-PD) coupling techniques in simulating large-scale geological hazards controlled by pre-existing fractures, this invention aims to solve the following specific technical problems:

[0010] 1. Existing coupling mechanisms are not well-suited to the "pre-existing fracture-settlement driven" scenario: Existing coupling methods (such as dynamic element transformation, sub-model iteration, etc.) are mostly designed for industrial material fracture or hydraulic fracturing scenarios. Their coupling logic does not fully consider the geological structural characteristics of the "pre-existing fracture" and the specific mechanical conditions of "pore water pressure load driven," making it difficult to accurately reproduce the fracture evolution details under complex loads while ensuring computational efficiency when simulating the entire process of large-scale ground subsidence accompanied by ground fissures.

[0011] 2. Simulation distortion caused by insufficient utilization of prior geological information in existing coupling architectures: Existing technologies mostly rely on stress-triggered passive conversion mechanisms, which cannot effectively integrate known prior fracture geometry and strike information in the initial stage of simulation. This makes it impossible for the model to accurately depict the "guiding role" of geologically weak structures on the initiation and propagation of ground fissures, resulting in deviations between simulation results and actual engineering geological conditions.

[0012] 3. Computational Bottlenecks and Numerical Oscillations Caused by High-Frequency Data Interaction: In large-scale simulations, existing coupling techniques involve cumbersome interface data interpolation, real-time mesh topology reconstruction, and stiffness matrix reorganization, resulting in a large amount of redundant algebraic calculations that are not directly related to the crack evolution mechanism. Furthermore, because transient transitions easily induce numerical oscillations and energy inaccuracies in the physical field, they struggle to meet the dual requirements of stability and efficiency for long-term evolution simulations at the engineering scale.

[0013] 4. Mismatch between existing rupture criteria and soil shear failure mechanisms: Most mainstream coupling methods adopt the bond elongation (tensile) failure criterion based on classical PD. This criterion cannot characterize the stress characteristics of soil and rock masses under uneven settlement loads, which are dominated by shear failure (i.e., it does not conform to soil and rock strength criteria such as Mohr-Coulomb), thus limiting the physical realism of ground fissure simulations.

[0014] In summary, the present invention aims to provide an efficient and accurate coupled simulation framework specifically for the aforementioned geological problems, in order to resolve the contradiction between continuous and discontinuous methods, and between computational efficiency and simulation accuracy in such problems. Summary of the Invention

[0015] The purpose of this invention is to provide a method and system for simulating ground fissures under pre-existing fracture geological conditions using coupled FEM-PD, in order to solve the problems mentioned in the background art, such as the insufficient adaptability of the existing coupling mechanism of the existing finite element-near field dynamics coupling technology to the "pre-existing fracture-settlement driven" scenario, the simulation distortion caused by insufficient utilization of geological prior information by the existing coupling architecture, the computational bottleneck and numerical oscillation caused by high-frequency data interaction, and the mismatch between the existing fracture criteria and the shear failure mechanism of rock and soil.

[0016] To achieve the above objectives, the present invention employs the following technical solution: In its first aspect, this invention proposes a method for simulating ground fissures under pre-existing fault geological conditions coupled with FEM-PD, comprising the following steps: S1. Construct the finite element-peripheral dynamics model and initialize the parameters; Based on the pre-existing fracture data obtained in the study area, the computational domain of the geological body is statically divided into a finite element region, a periphery dynamics region, and a coupling transition region between the two; The finite element region is discretized using a finite element mesh, and the periphery dynamics region is discretized using periphery dynamics material points; In the coupling transition region, the finite element mesh size is consistent with the spacing of the periphery dynamics material points, and the external load of the finite element region is transferred to the periphery dynamics region; S2. Finite element domain solution: The pore water pressure change caused by groundwater over-extraction is applied as an external load to the finite element domain. By solving the finite element control equations, the global displacement field and stress field of the finite element domain are obtained. S3, Unidirectional transfer of displacement boundary: The global displacement field obtained by finite element calculation is used as the displacement boundary condition and applied unidirectionally to the boundary of the near-field dynamic region. S4. Near-field dynamic region damage calculation: Based on the boundary conditions of the near-field dynamic region and the pore water pressure, the motion equations of the finite element-near-field dynamic model are obtained, and the motion equations of the near-field dynamic region are acquired. The adaptive dynamic relaxation method is used, combined with the motion equations of the near-field dynamic region, to perform time-domain integration calculation. Through step-by-step iteration, starting from the initial state, the displacement and stress of the material points at different times are calculated sequentially to simulate the motion process of the material points. The Mohr-Coulomb failure criterion is used to determine the breakage of the bond between two material points until all the bonds between material points are traversed, and local damage calculation is performed.

[0017] Preferably, the finite element region, the near-field dynamics region, and the coupling transition region between the two in S1 are as follows: The finite element region is used to simulate large-scale continuous deformation regions; The near-field dynamic region includes the pre-existing fault zone and the surrounding stress concentration area, used to simulate the evolution of ground fissures; The coupling transition region is defined as the finite element mesh adjacent to the near-field dynamic region. The coupling transition region includes coupling elements and coupling bonds. The finite element elements containing near-field dynamic material points are used as coupling elements, and only forces are applied to the finite element nodes in the coupling elements. The virtual bonds between the finite element nodes and the material points are used as coupling bonds, and only forces are applied to the near-field dynamic material points in the coupling bonds.

[0018] Preferably, step S2 is as follows: Calculate the stiffness matrix and load vector of all mesh elements. The stiffness matrix is ​​used to characterize the rigidity of the object under external force, and the load vector is used to characterize the external load. Based on the stiffness matrix and load vector, construct the finite element control equations and solve the finite element control equations using the double conjugate gradient stabilization method to obtain the displacement of each element node.

[0019] Furthermore, the solution of the finite element governing equations is as follows: Calculate the stiffness matrix and load vector Then, the stiffness matrices of all mesh elements and the load vectors are combined into a global stiffness matrix. With load vector Specifically:

[0020]

[0021]

[0022]

[0023]

[0024]

[0025] in, It is the external load vector; The finite element governing equations are solved using the double conjugate gradient stabilization method, specifically as follows:

[0026] in, It is a quality matrix; It is an acceleration vector; It is the global stiffness matrix; It is the nodal displacement vector of the finite element method; It is the external load vector; It is the strain operator (strain-displacement matrix); It is the elasticity matrix of the material (for isotropic elastic materials, Dependent on Young's modulus Compared to Poisson ); It is a partial derivative operator; It is a matrix of shape functions.

[0027] Preferably, step S3 is as follows: Finite element nodal displacement vectors As a displacement boundary condition, it is applied unidirectionally to the boundary material points of the near-field dynamic region; the displacement of the material points in the coupled transition region is obtained by interpolation from the adjacent finite element nodes through shape functions.

[0028] Preferably, the equations of motion for the finite element-peripheral dynamics model in S4 are as follows: Finite element nodal displacement vector As a displacement boundary condition, the finite element-near-field dynamics model at time step Internal force vector Internal force vectors of the finite element region and the internal force vector of the near-field dynamic region The combination yields, where Obtained by nonlocal integration; Represented as:

[0029] Among them, nodes and These are finite element nodes and near-field dynamic material points, respectively. It is the first n The displacement at each time step.

[0030] Furthermore, the internal force vector of the finite element region The details are as follows: The internal force vectors of a finite element region can be obtained by processing a continuous and regular region using the finite element method. Specifically:

[0031]

[0032]

[0033] in, The stiffness matrix of the finite element region; The internal force vector of the near-field dynamic region The details are as follows: The internal force vectors of the peridynamic region are obtained by using a peridynamic model to handle discontinuous regions. The calculation equation is as follows:

[0034] in, and It is the first Each time step node and scalar state of force density The unit vector state along the direction of the inter-point bond after deformation; node These are near-field dynamic particles or finite element nodes; It is the volume of a point mass.

[0035] Preferably, the equations of motion for the near-field dynamic region in S4 are as follows:

[0036] In the formula, It is the density of the material. yes acceleration vector; It is a point of matter Volume; It is the volume force density of an external force acting on a point of matter; It is a point of matter The near field region, , It is the radius of the near-field region; and They are matter points and The force density vector state at a point represents the material point in the model. and The interaction forces between them.

[0037] Preferably, in step S4, the Mohr-Coulomb failure criterion is used to determine the bond breaking status between two material points, as follows: At each time step, for each inter-material bond, the stress on its action surface is calculated; when the stress state satisfies the Coulomb failure criterion, that is, when the shear stress reaches the shear strength associated with the normal stress, it is determined that the inter-material bond has broken, and the local damage value is updated.

[0038] In a second aspect, this invention proposes a ground fissure simulation system coupled with FEM-PD under pre-existing fault geological conditions, comprising: The model discretization and static partitioning module is used to divide the geological body computational domain into a finite element region, a near-field dynamic region, and a coupling transition region between the two based on pre-existing fracture data. The parameter initialization module is used to input the required physical and computational parameters. This module injects the physical parameters required for the system calculation, including but not limited to geomechanical parameters such as soil density, Young's modulus, Poisson's ratio, cohesion, and internal friction angle, as well as numerical calculation parameters such as neighborhood range and time step.

[0039] The finite element solver is used to perform calculations in the finite element region. For the input pore water pressure variation load, it constructs and solves the global stiffness matrix equation and outputs the displacement field and stress field of the finite element region. The coupled data transfer module is used to apply the displacement field obtained by finite element calculation as displacement boundary condition to the boundary of the near-field dynamic region in one direction. The near-field dynamics solver is used to perform calculations in the near-field dynamics region. It performs nonlocal integral calculations based on the boundary conditions provided by the coupled data transfer module and uses the stress-based Coulomb failure criterion to determine the breakage of inter-point bonds in the material to simulate the evolution of ground fissures.

[0040] The model discretization and static partitioning module forms the topological foundation of the entire system, defining the scope of operation for the remaining modules. The parameter initialization module provides unified physical and numerical parameters for all subsequent calculation modules. In terms of the computational flow, the finite element solver performs calculations first under the drive of pore water pressure loads. Its output displacement field results are transmitted to the near-field dynamics solver via the coupled data transfer module, serving as input conditions for its calculations. Finally, the system integrates the results from the two solvers, outputting the full-field evolution process from continuous ground subsidence to discontinuous ground fissures.

[0041] By integrating the above processes, a complete, integrated simulation of the entire process, from large-scale ground subsidence (FEM simulation) to local ground fissure evolution (PD simulation), is achieved.

[0042] Compared with the prior art, the beneficial effects of the present invention are: (1) The method in this invention achieves a pioneering application: For the first time, the FEM-PD coupled model is applied to the simulation of ground subsidence-ground fissure evolution caused by groundwater extraction under pre-existing fracture conditions, which solves the problem that continuous and discontinuous methods in this field cannot be taken into account at the same time.

[0043] (2) The method in this invention achieves a balance between accuracy and efficiency: by using static partition coupling and unidirectional boundary transfer mechanism, while ensuring the accuracy of ground fissure simulation (see the embodiment, which is highly consistent with the results of the pure PD model), the computational efficiency is improved by about 75% compared with the pure PD model (the computation time is only 25.8%), making large-scale engineering simulation feasible.

[0044] (3) The method in this invention achieves a more accurate physical mechanism: the Coulomb failure criterion is adopted, which is more in line with the shear failure mechanism of rock and soil materials, and is more professional and accurate in ground fissure simulation than the general bond elongation criterion. Attached Figure Description

[0045] Figure 1 This is a flowchart of the ground fissure simulation method under pre-existing fault geological conditions coupled with FEM-PD in this invention; Figure 2 This is a schematic diagram of the finite element-peripheral dynamics coupling model in this invention; Figure 3 This is a schematic diagram of the geological body model in this invention; Figure 4 This is a schematic diagram of the conventional near-field dynamics model in this invention. Detailed Implementation

[0046] The technical solutions of 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.

[0047] Example 1: See Figure 1 A method for simulating ground fissures under pre-existing fault geological conditions using coupled FEM-PD is proposed, based on a finite element-near-field dynamics numerical simulation framework of "static partitioning-load-driven-unidirectional coupling". Its core working process involves pre-dividing the computational domain into a finite element region, a near-field dynamics region, and a coupling transition region between the two for a known pre-existing fault structure. This fixes the responsibilities of efficient computation (FEM) and accurate rupture (PD) from the outset, avoiding the computational overhead and uncertainties of dynamic conversion. The specific working process is as follows: Step 1: Model initialization and static partitioning.

[0048] Based on geological survey data, the locations of pre-existing faults are determined. Using this data, the computational domain of the entire geological body is statically and uniformly divided into three logical regions. The region dominated by large-scale, continuous deformation is designated as the finite element region, while fault zones and potential fracture areas are designated as the near-field dynamic region, with coupling transition zones established between them. Simultaneously, all geomechanical parameters are initialized.

[0049] (1) Model discretization: Construct a finite element-peripheral dynamics model. The structure of the finite element-peripheral dynamics model is as follows: Figure 2As shown, the geological model is discretized into a series of elements carrying physical information, and each element is divided into a finite element (FEM) region, a peri-field dynamics (PD) region, and a coupling transition region. Specifically: 1) Finite element region: The part used to simulate large-scale continuous deformation (i.e., ground settlement) is discretized using finite element mesh.

[0050] 2) Near-field dynamic region: This region is used to accurately simulate the initiation and propagation of ground fissures and is discretized using near-field dynamic material points. This region typically covers pre-existing fault zones and surrounding stress concentration areas.

[0051] 3) Coupling Transition Region: Located at the boundary between the two regions mentioned above, this region serves as a bridge for data exchange and displacement coordination. Finite element nodes and near-field dynamic material points coexist within this region. The core of the splicing coupling method lies in setting the finite element mesh adjacent to the near-field dynamic region as the coupling transition region, which is used as the near-field domain for the near-field dynamic material points. In the coupling transition region adjacent to the near-field dynamic region, the finite element mesh size maintains the same spacing as the near-field dynamic material points. Within the coupling transition region, finite element elements are divided into two categories based on whether they contain near-field dynamic material points: "regular elements" (without near-field dynamic material points) and "coupled elements" (containing near-field dynamic material points). Correspondingly, the "bonds" adjacent to the finite element region are divided into two categories based on whether they contain finite element nodes: "regular bonds" (without finite element nodes) and "coupling bonds" (containing finite element nodes) connecting two near-field dynamic material points. The key to achieving coupling between the near-field dynamics region and the finite element region lies in correctly defining the working mechanism of the "coupled element" and the "coupled bond": in the "coupled element", only the finite element nodes are subjected to force, while in the "coupled bond", only the near-field dynamic particles are subjected to force.

[0052] (2) Parameter initialization: Initialize all geomechanical parameters, including but not limited to: soil density, Young's modulus, Poisson's ratio, cohesion, internal friction angle, and other geomechanical parameters, as well as numerical calculation parameters such as neighborhood range and time step.

[0053] Step 2: Solve the finite element domain.

[0054] The calculations are performed in the finite element region using a finite element solver. It receives external pore water pressure variation loads and outputs the displacement and stress fields of the region by constructing and solving the global stiffness matrix equations.

[0055] (1) Calculation of global stiffness matrix: The basic process of the finite element method includes: structural discretization, element analysis, and global analysis.

[0056] First, as shown in formulas (1)-(6), calculate the stiffness matrix of each grid element in the discretized geological body model. Used to describe the rigidity of a material under external forces, and a load vector used to transmit external loads (such as gravity, mechanical loads, pore water pressure, etc.). Then, the stiffness matrices of all mesh elements and the load vectors are combined into a global stiffness matrix. With load vector .

[0057]

[0058] in, It is the external load vector; It is the strain operator (strain-displacement matrix); It is the elasticity matrix of the material (for isotropic elastic materials, Dependent on Young's modulus Compared to Poisson ); It is a partial derivative operator; It is a matrix of shape functions.

[0059] The pore water pressure changes caused by groundwater over-extraction are treated as external loads and directly applied to the finite element region. By solving the finite element governing equations, the global displacement and stress fields of the study area are obtained. This process clarifies that the driving source of the system is the effective stress change under seepage, rather than the complex flow of fluid within the fractures.

[0060] 2) Solve the finite element equations: The pore water pressure change is applied as an external load within the finite element region. The double conjugate gradient stabilization method is used to solve the control equation of the finite element model, as shown in formula (7), to obtain the displacement of each element node.

[0061] (7) In the formula, It is the mass matrix, which reflects the mass distribution characteristics of the structure and plays a key role in the structure's motion response in scenarios such as dynamic analysis. It is an acceleration vector that describes the acceleration changes at each node of the structure; It is the overall stiffness matrix, which reflects the structure's ability to resist deformation; It is a nodal displacement vector, representing the displacement state of each node in the structure; It is the external load vector.

[0062] The solution to this system of equations employs the double conjugate gradient stabilization method. Due to the coupled pore water pressure loads and complex boundary conditions, the assembled global stiffness matrix... KThese are typically large, sparse, asymmetric matrices. Traditional direct solution methods are computationally expensive for such matrices, while conventional iterative methods may encounter convergence difficulties. The biconjugate gradient stabilization method is an iterative algorithm specifically designed for solving asymmetric linear equation systems. By introducing a dual-projection mechanism and combining it with the smoothing technique of the generalized minimum residual method, it significantly improves the convergence stability and computational efficiency when solving complex problems like this model, thus ensuring the feasibility of large-scale simulations.

[0063] Step 3: Unidirectional transfer of displacement boundary.

[0064] Displacement field data calculated by a unidirectional finite element solver As a displacement boundary condition, it is precisely applied to the boundary of the near-field dynamic region. This invention uses the displacement field obtained from finite element method calculation as a displacement boundary condition and applies it unidirectionally to the near-field dynamic region; this efficient "FEM-guided PD" coupling method ensures global deformation compatibility while avoiding complex bidirectional iterative solutions, significantly improving computational efficiency.

[0065] Step 4: Near-field dynamic damage calculation and professional criteria.

[0066] The calculations are performed in the peri-field dynamics region using a peri-field dynamics solver. Based on boundary conditions, it performs nonlocal integral calculations and employs a stress-based Coulomb failure criterion to determine the breaking of inter-point bonds in the material, thereby simulating the evolution of ground fissures.

[0067] In the near-field dynamic region, nonlocal calculations are performed based on the transmitted displacement boundary conditions. The stress-based Coulomb failure criterion is used to determine the breakage of the "bonds" between material points. This criterion is more consistent with the shear failure mechanism of rock and soil, thus enabling a more professional and accurate simulation of the initiation and propagation path of ground fissures.

[0068] The core of the finite element-peripheral dynamics coupling process is the calculation of the internal forces of the coupled model. It is assumed that in the global system, the nodes... and If the finite element nodes and near-field dynamic particles are respectively, then the internal force vector of the coupled model can be expressed as shown in formula (8): (8) Combined with the stiffness matrix of the finite element region displacement vector of the element The finite element method is used to process continuous and regular regions to obtain the internal force vectors of the finite element region. As shown in formulas (9)-(11): (9) (10) (11) The internal force vectors of the peridynamic region are obtained by processing discontinuous regions (such as damaged or fractured areas) using a peridynamic model. The calculation equation is shown in formula (12): (12) In the formula, and It is the first Each time step node and The force density scalar state, where the nodes It can be a near-field dynamic particle or a finite element node.

[0069] Step 5: Integrated solution throughout the entire process.

[0070] An adaptive dynamic relaxation method is employed, combined with the motion equations of the coupled model, to perform time-domain integral calculations. Through iterative calculations, starting from the initial state, the displacement and stress of the material point at different times are obtained sequentially, thus describing the motion process of the material point. Finally, the stress-based Coulomb failure criterion is used to determine the fracture of the "bond" and calculate the local damage rate.

[0071] Experimental verification: To verify the applicability and numerical efficiency of the finite element-peripheral dynamics coupled model in simulating continuous ground subsidence and discontinuous ground fissures caused by groundwater over-extraction under pre-existing fracture conditions, this invention studies both continuous and discontinuous numerical simulation cases, and compares the results of the coupled model with those of the finite element model and the periphery dynamics model. The cases are mainly divided into four parts: setting up and initializing the geological body model, model construction and static partitioning, solution process and application of core formulas, and verification of model effectiveness and efficiency. Part 1: Geological Model Setup and Parameter Initialization; First, establish the computational model. For example... Figure 3 As shown, a three-dimensional geological model with dimensions of 600 m × 200 m × 200 m (length × height × thickness) was constructed. A pre-existing fault zone with a dip angle of 80° was included in the model. The model contains two aquifer systems: on the left side of the fault zone, the thicknesses of the upper and lower aquifers are 50 m and 150 m, respectively; on the right side, they are 150 m and 50 m, respectively. The simulated load condition is as follows: pore water pressure variation is applied only in the aquifer on the left side of the fault zone, with its value increasing linearly from 0 to 4.905 × 10⁻⁶ m. 5 Pa was used to simulate the effect of a 50 m drop in the groundwater level.

[0072] Based on engineering geological survey data, the geomechanical parameters are assigned to this model as shown in Table 1.

[0073] Table 1 Geomechanical parameters

[0074] Part Two: Model Building and Static Partitioning.

[0075] According to the core idea of ​​this invention, the above-mentioned geological bodies are statically partitioned and discretized: (1) Finite element region: located in the outer boundary region of the model. Triangular elements are used for discretization, and 120, 40 and 40 meshes are divided along the length, height and thickness of the model, respectively, with a mesh size of 5 m.

[0076] (2) Near-field dynamic region: covering the existing fracture zone and surrounding potential fracture zones. Specifically, it is set as a region with a length × height × thickness of 400 m × 100 m × 100 m in the model. It is discretized using material points with a layer spacing of 5 m, and the near-field range (neighborhood radius) is set to 3.015 times the layer spacing of the material points. To handle boundary effects, three layers of virtual boundary layer material points are set.

[0077] (3) Coupling transition region: This is the boundary between the finite element region and the near-field dynamics region, with a width of 30 m. Material points and finite element nodes in this region interact through "coupling bonds," serving as a bridge for data transfer.

[0078] Part Three: Solution Process and Application of Core Formulas.

[0079] The solution process of this invention is based on Figure 1 The calculation process is shown, and the calculation formulas listed in Example 1 are applied in detail. Furthermore, the solution for the near-field dynamic region employs the stress-based Coulomb failure criterion and the adaptive dynamic relaxation method; the core formulas and iterative process are explained in detail below.

[0080] (1) Finite element domain solution: The change in pore water pressure is used as the external load vector. The displacement field is applied to the finite element region. The global displacement field is solved using the following steps: 1) As shown in formulas (1) to (6), calculate the element stiffness matrix of each finite element. and load vector And assembled into a global stiffness matrix. K and global load vector F .

[0081] 2) Solve the finite element governing equations according to formula (7):

[0082] In the formula, It is a quality matrix. It is an acceleration vector. It is the global stiffness matrix. It is a nodal displacement vector. It is the external load vector.

[0083] (2) Unidirectional transmission of coupled data: The finite element nodal displacements obtained in step (1) As a displacement boundary condition, it is applied unidirectionally to the boundary material points of the near-field dynamic region. For material points within the coupling region, their displacements are obtained by interpolation from adjacent finite element nodes using shape functions.

[0084] (3) Solving the near-field dynamics region and determining damage: In the near-field dynamic region, nonlocal calculations are performed based on the displacement boundary conditions passed from step (2), specifically including: 1) Calculate the internal forces of the coupled system.

[0085] As shown in equations (8) to (12), the coupled model at time step Internal force vector Internal forces in the finite element region and internal forces in the near-field dynamic region Assembled. Among them, in formula (8) It can be calculated by nonlocal integration of formula (12).

[0086] Based on the aforementioned internal forces, such as Figure 4 As shown, the near-field dynamics equations of motion are described as follows:

[0087] In the formula, It is the density of the material. yes Acceleration vector. It is a point of matter The volume. It is the volume force density of an external force acting on a particle. It is a point of matter The near-field region is defined as: , It is the radius of the near-field region. and They are matter points and The force density vector state at a given point represents the material point in the deformed configuration. and The interaction forces between them. Matter points. and The original relative position and the deformed relative position are defined as the initial relative position vector state. and Deformation relative position vector state .

[0088] Force density vector state Defined as:

[0089] in, It is a scalar state of force density. For a small-strain linear elastic material, its expression is:

[0090] In the formula, and These are material parameters. It is the degree of volume expansion. It is a deflected extended state. m It is a weighted volume. It is an influence function used to control the effects within the near-field domain. Let the unit vector state along the direction of the "bond" after its deformation be defined as follows: ，

[0091]

[0092] =

[0093] )

[0094]

[0095]

[0096]

[0097]

[0098] Incorporating pore water pressure into the state-type near-field dynamics model:

[0099] In the formula, It is the fluid pressure coefficient. It is pore water pressure.

[0100] 2) Adaptive dynamic relaxation method: According to the adaptive dynamic relaxation method, the equations of motion of a matter point can be written as a series of ordinary differential equations by introducing virtual inertia and damping terms:

[0101]

[0102] Speed ​​updates:

[0103] Displacement update:

[0104] In the formula, For the first In the next iteration, since the displacement field at time t is unknown, the velocity formula cannot be used at the beginning of the iteration, but it can still be assumed that... and ,get:

[0105] Damping coefficient Dynamically updated:

[0106] The values ​​on the main diagonal of the local stiffness matrix are represented as:

[0107] in, It is a virtual diagonal density matrix (virtual mass matrix). It is the damping coefficient. It is the acceleration of a point mass. speed. It consists of near-field dynamic interaction forces and body forces. It is the first The resultant force vector at the time step. It is the first The speed of the step, It is the first The damping coefficient of the step, It is the time step. It is the inverse of the virtual mass matrix. It is the first The total force acting on a point mass. It is the first n The displacement at each time step.

[0108] 3) Damage assessment is performed using the Coulomb failure criterion.

[0109] At each time step, for each material point pair (bond), the stress on its action surface is calculated. When this stress state satisfies the Coulomb failure criterion (i.e., the shear stress reaches the shear strength associated with the normal stress), the bond is considered to have fractured, and the local damage value of the material is updated. For each material point pair (bond), the stress on its action surface is calculated as follows:

[0110]

[0111] Normal stress components:

[0112] Tangential stress components:

[0113] Shear strength:

[0114] The direction vector acting on the "bond" The traction force on the surface. The stress vector of a "bond" is defined as the average stress at the two material points connected by the "bond". For cohesion, It is the internal friction angle. When the modulus of tangential stress... or normal stress When the value is greater than 0, the bond is considered broken. Local damage value. Calculated based on the proportion of broken bonds:

[0115] in, =0 indicates that the key is complete. =1 indicates that the bond is broken. Conditional breakage for a single "bond":

[0116] (4) Integrated solution and output: The adaptive dynamic relaxation method is used to integrate the motion equations of the coupled system in the time domain, and the above steps are repeated until the simulation ends. Finally, the system outputs the complete process results from large-scale ground subsidence (characterized by the finite element displacement field) to local ground fissure evolution (characterized by the near-field dynamic damage field).

[0117] Part Four: Model Validity and Efficiency Testing.

[0118] (1) The comparison of continuous deformation results is shown in Table 2: the extreme value differences of the calculation results of the three models are small. Among them, the maximum difference in horizontal displacement is 0.09 m, accounting for 9.28% of the maximum horizontal displacement value; the maximum difference in vertical displacement is 0.03 m, accounting for 2.75% of the maximum vertical displacement value. The displacement extreme value of the coupled model is always within the extreme value range of the finite element model and the peri-field dynamics model, indicating that the simulation accuracy of the coupled model is better.

[0119] Table 2 Comparison of deformation extrema simulated by the finite element model, peridynamic model, and finite element-peridynamic coupling model.

[0120] (2) The comparison of discontinuous deformation results is shown in Table 3: the maximum horizontal displacement of the coupled model is consistent with that of the peri-field dynamic model, while the maximum vertical displacement is 0.01 m smaller than that obtained by the peri-field dynamic model. Similar to the aforementioned continuous case, the horizontal displacement of the coupled model shows a roughly smooth transition, demonstrating good numerical stability.

[0121] Table 3. Comparison of deformation extrema simulated by the peri-field dynamics model and the finite element-peri-field dynamics coupled model.

[0122] (3) Computational efficiency comparison is shown in Table 4: The computational efficiency of the finite element-peripheral dynamics coupled model and the periphery dynamics model are compared to evaluate the advantages of the coupled model in reducing computational complexity and saving computational costs. Both models are run on the same hardware platform, which is a macOS Sonoma operating system equipped with an 8-core CPU (4 performance cores + 4 efficiency cores), a 10-core GPU, and 24 GB of unified memory. Both models use the same discrete size, with a total of 241,875 material points and a time step of 400.

[0123] Table 4. Comparison of ground fissure characteristics simulated by the near-field dynamics model and the finite element-near-field dynamics coupled model.

[0124] The results show that, under the same hardware and discrete dimensions, the performance comparison between the coupled model of this invention and the pure peridynamic model is shown in Table 4. The finite element-peridynamic coupled model exhibits a significant computational efficiency advantage, with its total computation time being only 25.80% of that of the peridynamic model and its CPU usage being only 13.80% of that of the peridynamic model, resulting in a computational efficiency improvement of approximately 75%. This fully demonstrates the significant efficiency advantage of this invention in solving large-scale problems.

[0125] The above description is only for the purpose of helping to understand the method and core essence of the present invention, but the scope of protection of the present invention is not limited thereto. For those skilled in the art, any equivalent substitutions or modifications made to the technical solution and inventive concept disclosed in the present invention within the scope of the technology disclosed in the present invention should be covered within the scope of protection of the present invention. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A method for simulating ground fissures in pre-existing fractured geological conditions coupled with FEM-PD, characterized in that, Includes the following steps: S1. Construct the finite element-peripheral dynamics model and initialize the parameters; Based on the pre-existing fracture data obtained in the study area, the computational domain of the geological body is statically divided into the finite element region, the periphery dynamics region, and the coupling transition region between the two; The finite element region is discretized using a finite element mesh, while the near-field dynamics region is discretized using near-field dynamics material points. In the coupled transition region, the finite element mesh size is consistent with the near-field dynamic material point spacing, transferring external loads from the finite element region to the near-field dynamic region; S2, Finite element domain solution; The change in pore water pressure caused by groundwater over-extraction is applied as an external load to the finite element region. By solving the finite element control equations, the global displacement field and stress field of the finite element region are obtained. S3, Unidirectional transmission of displacement boundary; The global displacement field obtained by finite element calculation is used as the displacement boundary condition and applied unidirectionally to the boundary of the near-field dynamic region. S4. Damage calculation of the near-field dynamic region; Based on the boundary conditions of the near-field dynamic region and the pore water pressure, the motion equation of the finite element-near-field dynamic model is obtained, and the motion equation of the near-field dynamic region is obtained. An adaptive dynamic relaxation method is used, combined with the near-field dynamic region motion equation, to perform time-domain integral calculations. Through step-by-step iteration, starting from the initial state, the displacement and stress of the material points at different times are calculated sequentially to simulate the motion process of the material points. The Mohr-Coulomb failure criterion is used to determine the breakage of the bonds between two material points until all the bonds between material points are traversed, and local damage calculations are performed.

2. The method for simulating ground fissures under pre-existing fault geological conditions coupled with FEM-PD according to claim 1, characterized in that, The finite element region, the near-field dynamics region, and the coupling transition region between the two in S1 are specifically as follows: The finite element region is used to simulate large-scale continuous deformation regions; The near-field dynamic region includes the pre-existing fault zone and the surrounding stress concentration area, used to simulate the evolution of ground fissures; The coupling transition region is defined by setting the finite element mesh adjacent to the near-field dynamics region as the coupling transition region. The coupling transition region contains coupling elements and coupling bonds. Finite element elements containing near-field dynamic material points are used as coupling elements, and forces are applied only to the finite element nodes in the coupling elements. The virtual bond between the finite element node and the material point is used as the coupling bond, and the force is applied only to the near-field dynamic material point in the coupling bond.

3. The method for simulating ground fissures under pre-existing fault geological conditions coupled with FEM-PD according to claim 1, characterized in that, S2 is specifically as follows: Calculate the stiffness matrix and load vector of all mesh elements. The stiffness matrix is ​​used to characterize the rigidity of the object under external force, and the load vector is used to characterize the external load. Based on the stiffness matrix and load vector, construct the finite element control equations and solve the finite element control equations using the double conjugate gradient stabilization method to obtain the displacement of each element node.

4. The method for simulating ground fissures under pre-existing fault geological conditions coupled with FEM-PD according to claim 3, characterized in that, The solution to the finite element governing equations is as follows: Calculate the stiffness matrix and load vector Then, the stiffness matrices of all mesh elements and the load vectors are combined into a global stiffness matrix. With load vector Specifically: in, It is the external load vector; The finite element governing equations are solved using the double conjugate gradient stabilization method, specifically as follows: in, It is a quality matrix; It is an acceleration vector; It is the global stiffness matrix; It is the nodal displacement vector of the finite element method; It is the external load vector; It is a strain operator; It is the elasticity matrix of the material; It is a partial derivative operator; It is a matrix of shape functions.

5. The method for simulating ground fissures under pre-existing fault geological conditions coupled with FEM-PD according to claim 4, characterized in that, S3 is specifically as follows: Finite element nodal displacement vectors As a displacement boundary condition, it is applied unidirectionally to the boundary material points of the near-field dynamic region; the displacement of the material points in the coupled transition region is obtained by interpolation from the adjacent finite element nodes through shape functions.

6. The method for simulating ground fissures under pre-existing fault geological conditions coupled with FEM-PD according to claim 1, characterized in that, The equations of motion for the finite element-peripheral dynamics model in S4 are as follows: Finite element nodal displacement vector As a displacement boundary condition, the finite element-near-field dynamics model at time step Internal force vector Internal force vectors of the finite element region and the internal force vector of the near-field dynamic region The combination yields, where Obtained by nonlocal integration; Represented as: Among them, nodes and These are finite element nodes and near-field dynamic material points, respectively. It is the first n The displacement at each time step.

7. The method for simulating ground fissures under pre-existing fault geological conditions coupled with FEM-PD according to claim 6, characterized in that, The internal force vectors of a finite element region can be obtained by processing a continuous and regular region using the finite element method. ; The internal force vectors of the peridynamic region are obtained by using a peridynamic model to handle discontinuous regions. The calculation equation is as follows: in, and It is the first Each time step node and scalar state of force density The unit vector state along the direction of the inter-point bond after deformation; node These are near-field dynamic particles or finite element nodes; It is the volume of a point mass.

8. The method for simulating ground fissures under pre-existing fault geological conditions coupled with FEM-PD according to claim 6, characterized in that, The motion equations for the near-field dynamic region in S4 are as follows: In the formula, It is the density of the material. yes acceleration vector; It is a point of matter Volume; It is the volume force density of an external force acting on a point of matter; It is a point of matter The near field region, , It is the radius of the near-field region; and They are matter points and The force density vector state at a point represents the material point in the model. and The interaction forces between them.

9. The method for simulating ground fissures under pre-existing fault geological conditions coupled with FEM-PD according to claim 1, characterized in that, In S4, the Mohr-Coulomb violation criterion is used to determine the bond breaking status between two material points, as detailed below: At each time step, for each inter-material bond, the stress on its action surface is calculated; when the stress state satisfies the Coulomb failure criterion, that is, when the shear stress reaches the shear strength associated with the normal stress, it is determined that the inter-material bond has broken, and the local damage value is updated.

10. A ground fissure simulation system for pre-existing fault geological conditions coupled with FEM-PD applied to the method of any one of claims 1-9, characterized in that, include: The model discretization and static partitioning module is used to divide the geological body computational domain into a finite element region, a near-field dynamic region, and a coupling transition region between the two based on pre-existing fracture data. The parameter initialization module is used to input the required physical and calculation parameters; The finite element solver is used to perform calculations in the finite element region. For the input pore water pressure variation load, it constructs and solves the global stiffness matrix equation and outputs the displacement field and stress field of the finite element region. The coupled data transfer module is used to apply the displacement field obtained by finite element calculation as displacement boundary condition to the boundary of the near-field dynamic region in one direction. The near-field dynamics solver is used to perform calculations in the near-field dynamics region. It performs nonlocal integral calculations based on the boundary conditions provided by the coupled data transfer module and uses the stress-based Coulomb failure criterion to determine the breakage of inter-point bonds in the material to simulate the evolution of ground fissures.

Images