Efficient PD-FEM-FVM simulation and analysis method for construction rock mass stress-seepage coupling, and system

Abstract

Provided are an efficient PD-FEM-FVM simulation and analysis method for construction rock mass stress-seepage coupling, and a system. An embodiment comprises: separately establishing models of corresponding sizes for a solid layer and a fluid layer according to actual construction; applying initial boundary conditions on to the model as a whole; respectively performing calculation on a solid layer and a fluid layer on the basis of the applied initial boundary conditions, determining whether bonds of material points meet a destruction condition according to displacement circumstances of peridynamic material points in the solid layer, recording a local damage situation, and further describing fracturing of the rock mass; and obtaining the flow field distribution of the model as a whole according to the pressure of an FVM cell in the fluid layer, and thereby performing efficient construction-scale simulation of rock mass stress-seepage coupling.

Description

The present invention claims the priority of the Chinese patent application filed with the China National Intellectual Property7 Administration on November 30, 2022, with the application number 202211520123.9 and entitled "EFFICIENT SIMULATION AND ANALYSIS METHOD AND SYSTEM FOR STRESS-SEEPAGE COUPLING IN ENGINEERING ROCK MASS BASED ON PD-FEM-FVM", the entire content of which is incorporated herein by reference. TECHNICAL FIELD The present invention belongs to the technical field of numerical simulation of fluid-solid coupling in rock mass, and particularly relates to an efficient simulation and analysis method and system for stress-seepage coupling in engineering rock mass based on peridynamics (PD)-finite element method (FEM)-fmite volume method (FVM). BACKGROUND The statements in this section merely provide background technical information related to the present invention and do not necessarily constitute the prior art. In recent years, infrastructure construction such as railways, water conservancy and hydropower projects, and municipal projects has been booming, and China's transportation construction has ushered in unprecedented development opportunities. The increasing saturation of highway transportation has greatly promoted the construction of underground engineering represented by tunnels. Due to the presence of groundwater, geological disasters in tunnels and underground engineering are mostly problems of tire coupling of the stress field and the seepage field. Scientifically understanding the evolution mechanism of disasters is the theoretical basis and scientific basis for disaster prevention and control. With the rapid development of computer technology, numerical simulation has the advantages of strong visualization, good repeatability, and saving time and effort, and has become the main research means at present. The key to the simulation of stress-seepage coupling of engineering rock mass lies in the continuous-discontinuous evolution process. Peridynamics (PD) is one of the most advanced numerical calculation methods for solving discontinuous problems at present. It replaces the traditional differential model with an integral model of non-local action, and there is no problem of tip singularity. It is particularly suitable for simulating continuous-discontinuous problems. It is precisely because of the non-locality of PD that the calculation amount is relatively laige during the calculation, and the calculation efficiency is low, which is not conducive to large-scale calculations at the engineering scale. In solid calculation and fluid-solid coupling calculation, coupling PD with traditional continuum mechanics can give full play to the advantages of both, greatly improving the calculation efficiency, and is of great significance for realizing the calculation of fluid-solid coupling of rock mass at the engineering scale. Hie inventors have found that in the current coupling method of PD and FVM in fluid-solid coupling calculation, the fluid layer is calculated by adopting FVM, and the solid layer is all calculated by adopting PD. The coupling is realized by establishing a transition layer between the two to transfer information. Although this method improves the efficiency of fluid-solid coupling calculation to a certain extent, the solid layer is still all calculated by PD, the calculation efficiency of the solid portion is low, the amount of information interaction is large, and so the calculation efficiency can be further improved. SUMMARY In order to overcome the above deficiencies of the prior art, the present invention provides an efficient simulation and analysis method for stress-seepage coupling in engineering rock mass based on PD-FEM-FVM. The calculation efficiency can be further improved to achieve large-scale fluid-solid coupling calculations of rock mass at an engineering scale. In order to achieve the above objective, one or more embodiments of the present invention provide the following technical solutions: In a first aspect, an efficient simulation and analysis method for stress-seepage coupling in engineering rock mass based on PD-FEM-FVM is provided, including: establishing corresponding size models for a solid layer and a fluid layer separately according to actual engineering as rock mass calculation models; applying initial boundary conditions to the established rock mass calculation models; performing calculations for the solid layer and the fluid layer separately based on the applied initial boundary conditions, determining whether a bond of PD material points in the solid layer meets failure conditions according to displacement situations of material points, recording a local damage situation, and thereby describing a fracture situation of the rock mass; and obtaining a flow field distribution of an entire model according to pressure of an FVM control volume in the fluid layer, and thereby performing efficient simulation of the stress-seepage coupling in the rock mass at an engineering scale. As a further technical solution, an established model of a solid portion is divided into three regions, where a region where the rock mass is likely to fracture is a PD region, a region far away from a crack is a finite element region, and a portion connecting the two regions is a coupling region; the coupling region consists of connecting elements; and the connecting elements are finite element units embedded with the PD material points, the embedded material points are capable of being used to transfer force and displacement, and also serve as a virtual boundary layer of the PD region, and serve as objects for boundary condition application to weaken an influence of a boundary effect on a calculation result. As a further technical solution, a size of the FVM control volume in the fluid layer is die same as a size of a PD material point in the solid layer, and one finite element unit in the solid layer corresponds to a plurality of FVM control volumes. As a further technical solution, actual environmental conditions of the rock mass are equivalent to physical quantities as the initial boundary conditions, and the physical quantities include stress, displacement and water pressure. As a further technical solution, when performing the calculation for the solid layer: the PD region and the finite element region are connected by transferring the force and the displacement through the coupling region, and in coupling of PD and FEM for the solid portion, calculation of the force and the displacement is as follows: magnitudes of the transferred force and displacement are determined by using shape function interpolation according to relative positions of the embedded PD material points and nodes of the connecting elements. As a further technical solution, the fluid layer is calculated by adopting FVM, and in a coupling process between PD-FVM of the solid layer and the FVM of the fluid layer, a way of fluid pressure transfer is as follows: fluid pressure is transferred to the solid layer, at this time, a transition layer is controlled by the PD-FEM of the solid layer, and magnitudes of fluid pressure acting on the material points and the nodes of the finite element units are determined according to relative positions between the PD material points and the finite element units and nodes of a fluid grid. The solid layer is calculated by adopting an adaptive dynamic force relaxation algorithm to quickly balance initial in-situ stress. As a further technical solution, when performing the calculation for the fluid layer, a flow field in a saturated porous and fractured medium is described by using Darcy's law; and a pressure field of an entire region is obtained through a matrix operation by using the FVM. As a further technical solution, during the calculation for the fluid layer: a penneability coefficient of the solid layer is obtained, that is, solid deformation is transferred to a fluid grid, and at this time, a transition layer is controlled by FVM grid nodes of a fluid portion; and by determining PD material points and finite element units corresponding to control volumes of the grid nodes of the fluid portion, the solid deformation is transferred to the fluid grid by using permeability coefficients of these material points and nodes of the finite element units. Through the above calculations of the solid layer and the fluid layer, information such as the solid deformation and the fluid pressure can be obtained. In a second aspect, an efficient simulation and analysis system for stress-seepage coupling in engineering rock mass based on PD-FEM-FVM is provided, including: a model establishment module configured to: establish corresponding size models for a solid layer and a fluid layer separately according to actual engineering; a calculation module configured to: apply initial boundary conditions to the entire established models; perform calculations for the solid layer and the fluid layer separately based on the applied initial boundary' conditions, determining whether a bond of PD material points in the solid layer meets failure conditions according to displacement situations of material points, recording a local damage situation, and thereby describing a fracture situation of the rock mass; and a simulation module configured to: obtain a flow field distribution of an entire model according to pressure of an FVM control volume in the fluid layer, and thereby performing efficient simulation of the stress-seepage coupling in the rock mass at an engineering scale. The above one or more technical solutions have the following beneficial effects. (1) In one or more implementations of the present invention, by using the PD-FEM coupling method in the solid layer, the rock mass calculation model at the engineering scale is divided into the PD region, the finite element region, and the coupling region according to whether damage occurs. Compared with the pure PD method, the number of material points is greatly reduced, and the calculation workload is decreased. (2) In one or more implementations of the present invention, die fluid layer is calculated by adopting the FVM, which not only improves the calculation efficiency but also ensures the accuracy of the calculation results. (3) In one or more implementations of the present invention, the fluid-solid coupling calculation is realized by locally establishing the transition layer, which reduces the amount of information interaction between the solid layer and the fluid layer, further improves the efficiency of the fluid-solid coupling calculation, and thus enables efficient calculation of the fluid-solid coupling of the engineering rock mass. Hie additional advantages of the present invention will be partially given in the following description, partially become obvious from the following description, or be learned through the practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS The drawings of the specification that constitute a part of the present invention are used to provide a further understanding of the present invention. The schematic embodiments of the present invention and the descriptions thereof are used to explain the present invention and do not constitute an improper limitation to the present invention. FIG. 1 is a flowchart according to Embodiment 1 of the present invention. FIG. 2 is a PD-FEM hybrid model of a solid layer according to Embodiment 1 of the present invention. FIGS. 3(a)-3(b) are schematic diagrams of a PD-FEM-FVM fluid-solid coupling method according to Embodiment 1 of the present invention. FIG. 4 is a schematic diagram showing a corresponding relationship between an FEM element and a fluid grid according to Embodiment 1 of the present invention. FIG. 5 is a schematic diagram showing an FEM element bearing water pressure according to Embodiment 1 of the present invention. DETAILED DESCRIPTION It should be noted that the following detailed descriptions are all exemplary and are intended to provide further explanations of tire present invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by those of ordinary skill in the art to which the present invention belongs. It should be noted that the terms used herein are only for describing specific implementations and are not intended to limit the exemplary implementations according to the present invention. In the case of no conflict, the embodiments and the features in the embodiments of the present invention can be combined with each other. Embodiment 1 The present embodiment provides a numerical calculation method for stress-seepage coupling in engineering rock mass based on PD-FEM-FVM, wherein: first, a coupling model of PD and the FEM is established for a solid layer; second, a model is established for a fluid layer by adopting the FVM; and, a transition layer is established between a PD portion of the solid layer and a corresponding portion of the fluid layer, the transition layer stores a corresponding relationship between solid and fluid, thus realizing information transfer between fluid pressure and solid deformation. The method provided in the present embodiment includes the following steps. (1) Establishing a corresponding rock mass calculation model according to actual engineering, and dividing a solid layer into three regions, where a region where fracture may occur is a PD region, a region far away from a crack is a finite element region, and a portion connecting the two regions is a coupling region, which consists of connecting elements. The above connecting elements are finite element units embedded with the PD material points, the embedded material points are capable of being used to transfer force and displacement, and also serve as a virtual boundary layer of the PD region, and serve as objects for boundary condition application to weaken an influence of a boundary effect on a calculation result. A fluid layer consists of finite volume grids. A size of the finite volume grid in the fluid layer is the same as a size of the PD material point in the solid layer and smaller than that of FEM elements in the solid layer. Corresponding relationships between the PD material points in the solid layer, FEM nodes and corresponding finite volume grids in the fluid layer are recorded according to position coordinate relationships. During the recording, a program retrieves the corresponding relationships among the three according to position coordinates. Specifically, a PD-FEM hybrid model of the solid layer is shown in FIG. 2. PD-FEM-FVM fluid-solid coupling relationships are shown in FIGS. 3(a)-3(b), illustrating schematic diagrams of action of the fluid on the solid and action of the solid on the fluid, respectively. A schematic diagram of a corresponding relationship between an FEM element and a fluid grid is shown in FIG. 4. A schematic diagram of an FEM element bearing water pressure is shown in FIG. 5. The technical solutions of the present invention improve the efficiency of fluid-solid coupling calculation through the coupling calculation of PD and traditional continuum mechanic methods, so as to achieve large-scale calculations at an engineering scale. (2) Applying stress and displacement boundary conditions to the solid layer according to an actual engineering situation, and using an adaptive dynamic relaxation method for iterative calculation. In the solid layer, a region where an object is likely to fracture is a PD region, and a region far away from a crack is a finite element region. They are connected through interface elements, that is, force and displacement are mutually transmitted. Ilie interface elements are finite element units embedded with the PD material points, the embedded material points are capable of being used to transfer force and displacement, and also serve as a virtual boundary layer of the PD region, and serve as objects for boundary condition application to weaken an influence of a boundary effect on a calculation result. The number of the embedded material points is related to a size of a PD neighborhood range. Then, short-range repulsive force is introduced in PD. Therefore, an expression of coupling force is: m i=l Where, fcp is coupling force generated by all material points embedded in a certain interface element, m is the number of material points in the interface element, fl is bond force received by an embedded material point i, fr is short-range repulsive force received by the embedded material point i, specifically: r / n + E r cs >\ fi = j + min{°,y (1*7 + fl) ~ dvx' Where, ds = min{0.9|x — x'\, 1.35(rs + rs')} Where, rs is a radius of the material point. The coupling force is only distributed to finite element nodes on an interface between the PD region and the FEM region. Therefore, for two-dimensional (2D) problems, the coupling force distributed to the finite element nodes on the interface between the PD region and the FEM region is: C = ¥ / = 1,2,3,4 Where, Nt is a shape function of an element node I on the interface. Displacement of an embedded PD material point is obtained by interpolating displacement of the interface element node through the shape function: In order to achieve an efficient solution for static problems, by assembling a global matrix of a coupled system and using the adaptive dynamic relaxation method, the PD region and the FEM region are solved simultaneously within a unified framework, which speeds up convergence and greatly improves calculation efficiency. A discrete PD equation of motion can be written as: ft1 o o D-1 PS fS Where, Up is a displacement vector of a PD material point, Up and Up are second-order and first-order derivatives of displacement Up with respect to time respectively, i.e., acceleration and velocity vectors; a vector Fp is a sum of the internal force and external force received by a PD system; a superscript p represents a variable related to the PD region; single and double underlines represent variables located outside and inside an interface element respectively; cn represents a damping coefficient at an nth time step; and D is a virtual diagonal density matrix, and its diagonal elements are: Where is a stiffness matrix of the PD system. Similarly, a finite element equation can be expressed as: Where, a superscript f represents a variable related to the FEM region; and M is a diagonal mass matrix. The mass matrix can be approximated as: M = Im Where, I is a unit matrix. A mass vector m is: m = Am Where, A is an assembly operator. A component of m is: 7=1 Where, k^ is a component of a unit stiffness matrix. is an unbalanced force vector received by the finite element nodes at the nth time step and is expressed as: F^= (Ff 1 - (Ff 'I n \1extJn \ intJn Where, is an external force vector received by an FEM system at the ntfl time step; and «t)n is an internal force vector, and an internal force at any finite element node I is: . e J / W«),= 2 / fX) +fF - e J / node / is in a conventional FEM region node / is at an interface between PD and FEM regions Where, [■], represents a finite element unit related to a node I, is a stiffness matrix of an FEM element, is a displacement vector of the FEM element, and is coupling force acting on the finite element nodes located at the interface between the PD region and the FEM region. Finally, an equation set for the coupled system can be obtained as: U +cnU = M-1 F L—J__ I I—J I—J Where, D 0 0 0 M 0 OOM { / 7? f^ F^} In the adaptive dynamic relaxation method, a current optimal damping coefficient needs to be calculated at each time step to speed up the convergence. A damping coefficient c at the nth time step is: l(^n)T^Un Cn J (UnYUn Where, Kf is a diagonal local stiffness matrix, and its expression is: Fl — Fl In the present embodiment, the time step At is taken as 1. In order to obtain a steady-state solution, a central-difference explicit integration is adopted for calculation, and a recurrence formula is: (2 — cn) (7 1 + 2M~1F y =__________2__________ (2 + Cn) u =u +u , wn+l ^n+-| When t = 0, an iteration starts with: £1 2M Meanwhile, to consider the interaction between the solid deformation and the fluid pressure, according to the principle of effective stress, a finite element equation of motion is modified. The finite element equation of motion based on the principle of effective stress in a saturated porous medium is: U + cnU = - aP) Where, U is a displacement vector of the finite element nodes, U and U are second-order and first-order derivatives of displacement U with respect to time respectively, i.e., acceleration and velocity vectors; P is a water pressure vector borne by nodes of a finite element unit; cn represents a damping coefficient at the nth time step; a is a Biot coefficient; M is a diagonal mass matrix; and F is a force vector. A basic PD equation in the solid layer is modified. A relationship between a scalar force density state and a scalar effective force density state in PD is: {mx t total 3 3 m O)X itotal~^aP-- / 2D m Where, ttotai is the scalar force density state, is the scalar effective force density state, a is the Biot coefficient, p is a pore pressure, m is a weighted volume scalar function, and w is an influence function. Therefore, a 2D PD equation of motion based on the principle of effective stress in the saturated porous medium is: pii(x, t) = I (T[x, t](x' — x} — T[x', t](x — x'})dVx> JHx C wx wx — 2a I ip — H (x — x) — p —H(x — x Hc / Vv' + b(x, t) J tn ..................' 1 * T7X..................' f x v x JHX Where, p is a material density, b is a body force density, Hx is a neighborhood of a material point x, Y is a deformation, and T_ is a force vector state. Meanwhile, the PD portion of the solid layer is divided into three portions: a solid domain Qs, a transition domain Qt, and a fracture domain fy.Using a PD damage field as the basis for division, two indicators are set to identify the three regions, and a resulting linear index function is as follows: a — q c2 C1 C2 — u Xr ~ ~ c2 C1 Where, q and q are two division indicators, with q = 0.2 and q = 0.4. By interpolating permeabilities of the solid domain and the fracture domain using a linear index function, solid permeability7 of the transition domain can be obtained: / C = XrK + Xfy where kr and kf are the permeabilities of the solid domain and the fracture domain, respectively. The permeability of the fracture region is calculated using the cubic law: 1 2 — ..........................h ? 12M Where, b is a fracture aperture. In an iterative solution process, whether a bond of material points in the PD region meets failure conditions is determined, and a local damage situation is recorded. The failure conditions are a determination of integrity of the bond of the material points expressed by a critical elongation rate: (.[t.f] = (1 VS<S» ; t0 otherwise Where, s0 is an elongation rate of a given bond of the material points, s is an elongation rate of the bond of the material points, which is expressed as s = where, r / is relative displacement between any two material points, and f is a relative elongation rate between any two material points. That is, when tensile deformation s of the bond of the material points exceeds the given limit value s0, the bond breaks. At this time, there is no longer an interaction relationship between two interacting material points. Local damage is defined as a ratio of the number of remaining intact bonds to the initial number of bonds after bonds of the material points are broken, which is expressed as: ,.,, , f„ <p(i, t) = 1--; Where, is a volume of the material point i. It is noted that, 0 <<p <1, where, 0 represents an intact state, 1 represents a completely damaged state, and values between 0 and 1 are quantitative representations of a local damage degree. Based on the above calculations, a damage situation of the PD portion in the solid layer can be obtained, and permeability of the PD portion can be updated according to a breakage situation of the bonds between material points. Since no fracture occurs in the FEM region, permeability of an FVM control volume corresponding to the FEM region remains unchanged and is initial permeability. (3) Applying fluid pressure boundary conditions, and according to a corresponding relationship between the fluid grid and a solid grid, transferring permeability of the solid layer to nodes of the fluid grid, in a finite element portion of the solid layer, directly transferring the permeability to the corresponding FVM control volume in the fluid layer, and for the PD portion of the solid layer, transferring the permeability according to a relative position relationship between the PD material points and finite volume control volumes in the fluid layer. A flow field in a saturated porous and fractured medium is described by using Darcy's law, which can be expressed as: K u =--VP Where: is fluid viscosity, K is a permeability tensor of a porous medium, and P is fluid pressure. For a homogeneous isotropic body, K = kl. According to the relative position relationship, water pressure borne by the PD material points in a solid portion is equal to water pressure at corresponding nodes of the fluid grid. Each finite element unit in the solid portion corresponds to multiple nodes of the fluid grid. Water pressure within each unit is obtained by averaging water pressures at the nodes of the fluid grid within the unit, which is calculated by adopting the following formula: i = 1,2,3 Where, Pe is water pressure within a finite element unit in the solid portion, pt is the water pressure at the node of the fluid portion grid, and n is the number of fluid grids contained within the finite element unit. For a 2D finite element unit with a certain thickness, water load home on each node within the unit can be calculated by using the following formula: Pe ■ Ih Pt = —----, i = 1,2,3,4 4 Where, Pt is water load borne by a unit node, I is a unit width, and h is a unit thickness. Moreover, according to positions of four nodes within the unit, a sign of water load in x and y directions on each node can be determined. That is, signs of water loads in the x and y directions on nodes 1, 2, 3, and 4 are negative-negative, positive-negative, positive-positive, and negative-positive in sequence. For element nodes shared by multiple units, water load thereon is obtained by superposition. When the solid deformation and failure are transferred to the fluid grid, it is mainly reflected in change and transfer of solid permeability. At this time, the transition layer is controlled by the nodes of the fluid grid. Similar to transfer of fluid pressure to the solid, by determining PD material points and finite element nodes corresponding to the nodes of the fluid grid, the solid deformation is transferred to the fluid grid by using permeability coefficients of the material points and finite element nodes. Through the above calculations, the characteristics of the solid deformation, damage, and the fluid flow field can be obtained. (4) Determining whether a specified time step has been reached. If the specified time step is reached, the calculation is ended; and if not, the process proceeds to the next iterative step for calculation. The detailed calculation process of the technical solution in the present embodiment is shown in FIG. 1, which specifically includes: 1) performing initialization setting on basic variables and corresponding parameters; 2) determining geometric dimensions and material parameters of a corresponding numerical calculation model according to a simulation object, and dividing a PD region and an FEM region according to an expected failure mode, determining the size of a finite element grid, the size of a PD material point, and the range of a PD neighborhood, discretizing each region and embedding the corresponding number of material points within interface elements, and numbering discrete finite element nodes and PD material points, and determining node composition of each finite element unit; 3) determining the total number of calculation time steps and a time step size, where since an adaptive dynamic relaxation method is adopted as a solution strategy, a time step is generally taken as 1; 4) determine corresponding relationships between solid material points and finite element nodes and nodes of a fluid grid according to relative position relationships; 5) first performing a calculation of a solid portion, determining other material points within a neighborhood range of each PD material point, that is, determining group members of each material point, and storing this information; 6) calculating each material point for surface correction, that is, calculating a surface correction factor; 7) calculating an overall stiffness matrix and a mass matrix of the finite elements, and a stable mass matrix of PD, and assembling to form a global mass matrix of a coupled system; 8) applying initial boundary conditions and load conditions to a discrete model, and calculating an external force vector; 9) displacement transfer: interpolating displacement of finite element nodes on interface elements at a previous time step to obtain displacement of embedded material points within the interface elements, which are used as boundary conditions of a PD region, where displacement of finite element nodes on interface elements at an initial time step are set as initial values, i.e. zero; 10) calculating pairwise force and short-range repulsive force of each material point according to the boundary conditions of the PD region; 11) force transfer: summing up the pairwise force and the short-range repulsive force received by each material point within the interface elements, i.e., coupling force, and distributing it to finite element nodes on corresponding interface elements located at an interface between two regions; 12) calculating unbalanced force in a finite element region, and combining the pairwise force and the short-range repulsive force in the PD region to assemble and form an overall force matrix of the coupled system; 13) solving for displacement of the coupled system according to the adaptive dynamic relaxation method; 14) calculating an elongation rate of a bond between the PD material points; and if it exceeds a critical elongation rate, determining that the bond is broken, and if it does not exceed the critical elongation rate, determining that the bond is not broken; 15) determining whether all the material points have been traversed in the above steps; and if all the material points are traversed in the above steps, proceeding to the next time step for calculation, and otherwise, returning to step 14) to continue the calculation; 16) calculating the permeabilities of the finite element nodes and the PD material points at this time; 17) transferring the permeability of a solid layer to the fluid grid according to the corresponding relationship between the solid layer and the fluid grid established in step 4); 18) calculating a matrix coefficient corresponding to each fluid grid node separately; 19) assembling an overall calculation matrix of fluid and solving for a flow field distribution of the fluid; and 20) determining whether the calculation has reached a specified time step; and if the calculation reaches the specified time step, the calculation is ended and the program is exited, and if the calculation does not reach the specified time step, the process returns to step 8) to continue the calculation. Embodiment 2 An objective of the present embodiment is to provide a computer device, including a memory, a processor, and a computer program stored on the memory and executable on the processor, where when the processor, when executing the program, implements steps of the method described above. Embodiment 3 An objective of the present embodiment is to provide a computer-readable storage medium. A computer-readable storage medium having a computer program stored thereon, where the program, when executed by a processor, performs steps of the method described above. Embodiment 4 An objective of the present embodiment is to provide an efficient simulation and analysis system for stress-seepage coupling in engineering rock mass based on PD-FEM-FVM, including: a model establishment module configured to: establish corresponding size models for a solid layer and a fluid layer separately according to actual engineering; a calculation module configured to: apply initial boundary conditions to the entire established models; perform calculations for the solid layer and the fluid layer separately based on the applied initial boundaiy conditions, determining whether a bond of PD material points in the solid layer meets failure conditions according to displacement situations of material points, recording a local damage situation, and thereby describing a fracture situation of the rock mass; and a simulation module configured to: obtain a flow field distribution of an entire model according to pressure of an FVM control volume in the fluid layer, and thereby performing efficient simulation of the stress-seepage coupling in the rock mass at an engineering scale. The above solid layer portion: it includes portions such as model establishment, region division, boundary condition application, and calculation analysis. A calculation model of a corresponding size is established according to actual engineering, and divided into a PD region, a finite element region, and a coupling region. Stress and displacement boundary conditions are applied according to an actual situation, and calculation is performed by using an adaptive dynamic relaxation method to obtain stress, displacement, and damage situations; and the fluid layer portion: it includes portions such as model establishment, region division, boundary condition application, and calculation analysis. A flow field in a saturated porous and fractured medium is described by using Darcy's law, and a pressure field of an entire region is obtained through a matrix operation by using a FVM The transition layer portion is responsible for mutual transfer between solid deformation and fluid pressure. When the solid deformation is transferred to the fluid grid, the transition layer is controlled by the fluid portion. By determining PD material points and FEM nodes corresponding to FVM control volumes of the fluid portion, the solid deformation is transferred to the fluid grid by using the permeability of the material points and the FEM nodes. When the fluid pressure is transferred to the solid layer, the transition layer is controlled by the solid layer (PD-FEM). Similarly, a magnitude of fluid pressure received by the solid portion is determined according to relative positions among the PD material points, FEM nodes, and the fluid FVM control volumes. Calculation results (stress, displacement, damage, and pressure) of each iterative step are stored, and the calculation results are visually displayed in the form of a contour plot. The steps involved in the devices of Embodiments 2, 3, and 4 above correspond to those in the method of Embodiment 1. For specific implementations, reference can be made to the relevant description part of Embodiment 1. The term "computer-readable storage medium" should be understood as including a single medium or multiple media that contain one or more instruction sets. It should also be understood as including any medium that is capable of storing, encoding, or carrying an instruction set for execution by a processor and enabling the processor to execute any method of the present invention. Those skilled in the art should understand that the modules or steps of the present invention described above can be implemented by a general-purpose computer device. Optionally, they can be implemented by program codes executable by a computing device. Thus, they can be stored in a storage device and executed by the computing device, or they can be separately fabricated into individual integrated circuit modules, or multiple modules or steps among them can be fabricated into a single integrated circuit module for implementation. The present invention is not limited to any specific combination of hardware and software. Although the specific implementations of the present invention have been described in combination with the drawings above, it is not a limitation on the scope of protection of the present invention. Those skilled in the art should understand that, based on the technical solutions of the present invention, various modifications or variations that can be made by those skilled in the art without creative efforts still fall within the scope of protection of the present invention.

Claims

1. An efficient simulation and analysis method for stress-seepage coupling in engineering rock mass based on peridynamics (PD)-finite element method (FEM)-finite volume method (FVM), comprising:establishing corresponding size models for a solid layer and a fluid layer separately according to actual engineering as rock mass calculation models;applying initial boundary conditions to the established rock mass calculation models;performing calculations for the solid layer and the fluid layer separately based on the applied initial boundary conditions, determining whether a bond of PD material points in the solid layer meets failure conditions according to displacement situations of material points, recording a local damage situation, and thereby describing a fracture situation of the rock mass; andobtaining a flow field distribution of an entire model according to pressure of an FVM control volume in the fluid layer, and thereby performing efficient simulation of the stress-seepage coupling in the rock mass at an engineering scale.

2. Tire efficient simulation and analysis method for stress-seepage coupling in engineering rock mass based on PD-FEM-FVM according to claim 1, wherein an established model of a solid portion is divided into three regions, wherein a region where the rock mass is likely to fracture is a PD region, a region far away from a crack is a finite element region, and a portion connecting the two regions is a coupling region;the coupling region consists of connecting elements; andthe connecting elements are finite element units embedded with the PD material points, the embedded material points are capable of being used to transfer force and displacement, and also serve as a virtual boundary layer of the PD region, and serve as objects for boundary condition application to weaken an influence of a boundary effect on a calculation result.

3. The efficient simulation and analysis method for stress-seepage coupling in engineering rock mass based on PD-FEM-FVM according to claim 2, wherein a size of the FVM control volume in the fluid layer is the same as a size of a PD material point in the solid layer, and one finite element unit in the solid layer corresponds to a plurality of FVM control volumes.

4. The efficient simulation and analysis method for stress-seepage coupling in engineering rock mass based on PD-FEM-FVM according to claim 1, wherein actual environmental conditions of therock mass are equivalent to physical quantities as the initial boundary conditions, and the physical quantities comprise stress, displacement and water pressure.

5. The efficient simulation and analysis method for stress-seepage coupling in engineering rock mass based on PD-FEM-FVM according to claim 2, wherein when performing the calculation for the solid layer:the PD region and the finite element region are connected by transferring the force and the displacement through the coupling region, and magnitudes of the transferred force and displacement are determined by using shape function interpolation according to relative positions of tine embedded PD material points and nodes of the connecting elements.

6. The efficient simulation and analysis method for stress-seepage coupling in engineering rock mass based on PD-FEM-FVM according to claim 5, wherein in the calculation for the solid layer, when the transferred force is determined, fluid pressure is transferred to the solid layer, at this time, a transition layer is controlled by PD-FEM of the solid layer, and magnitudes of fluid pressure acting on the material points and the nodes of the finite element units are determined according to relative positions between the PD material points and the finite element units and nodes of a fluid grid.

7. The efficient simulation and analysis method for stress-seepage coupling in engineering rock mass based on PD-FEM-FVM according to claim 1, wherein when performing the calculation for the fluid layer, a flow field in a saturated porous and fractured medium is described by using Darcy's law; anda pressure field of an entire region is obtained through a matrix operation by using a FVM; and preferably, during the calculation for the fluid layer:a permeability coefficient of the solid layer is obtained, that is, solid deformation is transferred to a fluid grid, and at this time, a transition layer is controlled by FVM grid nodes of a fluid portion; andby determining PD material points and finite element units corresponding to control volumes of the grid nodes of the fluid portion, the solid deformation is transferred to the fluid grid by using permeability coefficients of these material points and nodes of the finite element units.

8. An efficient simulation and analysis system for stress-seepage coupling in engineering rock mass based on peridynamics (PD)-fmite element method (FEM)-fmite volume method (FVM),comprising:a model establishment module configured to: establish corresponding size models for a solid layer and a fluid layer separately according to actual engineering;a calculation module configured to: apply initial boundary conditions to the entire established models; andperform calculations for the solid layer and the fluid layer separately based on the applied initial boundary conditions, determine whether a bond of PD material points in the solid layer meets failure conditions according to displacement situations of material points, record a local damage situation, and thereby describe a fracture situation of the rock mass; anda simulation module configured to: obtain a flow field distribution of an entire model according to pressure of an FVM control volume in the fluid layer, and thereby performing efficient simulation of the stress-seepage coupling in the rock mass at an engineering scale.

9. A computer device, comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, wherein when the processor, when executing the program, implements steps of the method according to any one of claims 1-7.

10. A computer-readable storage medium having a computer program stored thereon, wherein the program, when executed by a processor, performs steps of the method according to any one of claims 1-7.INTERNATIONAL SEARCH REPORT International application No. PCT / CN2023 / 095^9A. CLASSIFICATION OF SUBJECT MATTER G06F 30 / 28(2020.01 )i According to International Patent Classification (IPC) or to both national classification and IPC B. FIELDS SEARCHED Minimum documentation searched (classification system followed by classification symbols) IPC: G06F Documentation searched other than minimum documentation to the extent that such documents are included in the fields searched Electronic data base consulted during the international search (name of data base and, where practicable, search terms used) VEN, CNABS, CNTXT, WOTXT, USTXT, EPTXT, CNKI, IEEE: g#, S±, PD, FEM, FVM, frWllg, WIWx, W , SBffl'H, IS pHt©, Ml#, Hff, 54 RK, peridynamics, fluid-solid coupling, coupling, finite volume, rock, fluid, solid C. DOCUMENTS CONSIDERED TO BE RELEV ANT Category* Citation of document, with indication, where appropriate, of the relevant passages Relevant to claim No. Y CN 113758848 A (SHANDONG UNIVERSITY) 07 December 2021 (2021-12-07) description, paragraphs 0002-0041, 0050-0087, and 0090-0110, and figures la and lb 1-10 Y A LIU, Wenyang et al. "A Coupling Approach of Discretized Peridynamics With Finite Element Method" Comput. Methods Appt. Meeh. Engrg., 31 December 2012 (2012-12-31), text, pages 166-167, and abstract CN 112131802 A (SHANDONG UNIVERSITY) 25 December 2020 (2020-12-25) entire document 1-10 1-10 A CN 113761760 A (SHANDONG UNIVERSITY) 07 December 2021 (2021-12-07) entire document 1-10 A (YU, Yangtian et al.). (1 lybrid Model of Peridynamics and Finite Element Method Under Implicit Schemes)" (Journal of Zhejiang UniversityfEngineering Science)), No. 07, 31 December 2017 (2017-12-31), entire document 1-10 | / 1 Further documents are listed in the continuation of Box C. | S | See patent family annex. * Special categories of cited documents: “T” later document published after the international filing date or priority “A” document defining the general state of the art which is not considered date and not in conflict with the application but cited to understand the to be of particular relevance principle or theory underlying the invention “D” document cited by the applicant in the international application “X” document of particular relevance; the claimed invention cannot be -E” earlier application or patent hut published on or after the international considered novel or cannot he considered to involve an inventive step filing date when the document is taken alone ■‘L” document which may throw doubts on priority claim(s) or which is "Y” document of particular relevance; the claimed invention cannot be cited to establish the publication date of another citation or other considered to involve an inventive step when the document is special reason (as specified) combined with one or more other such documents, such combination “O” document referring to an oral disclosure, use, exhibition or other being obvious to a pet son skilled in the art means document member of the same patent family “P" document published prior co the international filing date but later than the priority date claimed Date of the actual completion of the international search 08 September 2023 Date of mailing of the international search report 13 September 2023 Name and mailing address of the IS A / CN China National Intellectual Property Administration (ISA / CN) China No. 6, Xitucheng Road, Jimenqiao, Haidian District, Beijing 100088 Authorized officer Telephone No.Form PCT / ISA / 210 (second sheet) (July 2022)TRANSLATION INTERNATIONAL SEARCH REPORT International application No. PCT / CN2023 / 095869C. DOCUMENTS CONSIDERED TO BE RELEVANTCategory* Citation of document, with indication, where appropriate, of the relevant passages Relevant to claim No. A )• i*f ViZHANG. Qing et al.). (Hybrid Modeling Methods of Pendynamics and Finite Element Method)" (Chinese Journal of Computational Mechanics), No. 04,15 August 2016 (2016-08-15), entire document 1-10 A RichardW. Macek el al. "Peridynamics Via Finite Element Analysis" Available noline at www.sciencedirect.com, 31 December 2007 (2007-12-31), entire document 1-10 A US 2006129366 Al (SHAW, Gareth) 15 June 2006 (2006-06-15) entire document 1-10 — —INTERNATIONAL SEARCH REPORT Information on patent family members International application No. PCT / CN2023 / 095869Patent document cited in search report Publication date (day / month / year) Patent family member(s) Publication date (day / month / year) CN 113758848 A 07 December 2021 None CN 112131802 A 25 December 2020 None CN 113761760 A 07 December 2021 None US 2006129366 Al 15 June 2006 EP 1825411 A2 29 August 2007 EP 1825411 Bl 19 September 2018 NO 20073486 L 14 September 2007 NO 341648 Bl 18 December2017 WO 2006064379 A2 22 June 2006 WO 2006064379 A3 21 September 2006 CA 2596535 Al 22 June 2006 CA 2596535 C 25 March 2014 EA 200701279 Al 28 February 2008 EA 010^8 Bl 30 December 2008 US 7707018 B2 27 April 2010

Images