A numerical simulation method and apparatus for surface deformation evolution based on near-field dynamics

Abstract

This invention provides a numerical simulation method for surface deformation evolution based on near-field dynamics, comprising: discretizing the geological body model and virtual boundary layer to be calculated into multiple material points to obtain an initial geological body model; obtaining the neighborhood of each material point in the initial geological body model; obtaining the internal force experienced by the material point at the current time step based on the pore water pressure change at the current time step, as the stress on the material point at the current time step; obtaining the displacement of the material point at the current time step using an adaptive dynamic relaxation method based on the stress on the material point at the current time step; obtaining the stress on the material point at the current time step based on the displacement on the material point at the current time step; and determining the fracture status of the virtual bond using the Mohr-Coulomb failure criterion based on the stress on the material point at the current time step, thereby obtaining the local damage value of the material point at the current time step. This invention also provides another scheme capable of revealing the physical mechanism of ground subsidence disasters transforming into ground fissures.

Description

[0001] Technical field

[0002] The present application relates to the geological disaster field, in particular to a surface deformation evolution numerical simulation method and device based on near-field dynamics. BACKGROUND

[0003] Land subsidence is a geological phenomenon of loss of ground elevation. Under special geological structure conditions, differential land subsidence will evolve into ground fissures, which will cause serious threats to existing infrastructure, linear rail transit and underground space development and utilization, and restrict the sustainable development of economy and society. The evolution and development of land subsidence-ground fissures are affected by a variety of factors, and the differences in geological structure, soil lithology and groundwater exploitation cause different evolution characteristics and mechanisms of land subsidence. Understanding the evolution of land subsidence and the initiation and expansion mechanism of ground fissures is of great significance to risk management and sustainable development of underground resources.

[0004] Numerical simulation method is an important means to study the evolution and development process of land subsidence-ground fissures. The existing ground fissure simulation methods are mainly based on the finite element method and the finite difference method of continuous medium mechanics theory. However, the process of land subsidence disaster into ground fissure is a nonlinear and discontinuous time-space problem, and the above methods have contradictions between continuous theory and discontinuous reality. In addition, when using non-continuous medium mechanics theory, such as discrete element method, interface element method and extended finite element method, on the one hand, it is necessary to know in advance whether the crack exists and its position and size, and on the other hand, after the crack expands, the grid needs to be re-divided or additional functions need to be introduced to describe the discontinuous problem, which is difficult to realize the automatic and accurate simulation of complex cracks.

[0005] Near-field dynamics is based on the idea of non-local action, and describes the mechanical behavior of matter by solving spatial integral equations applicable to discontinuous bodies. The near-field dynamics method has a unified expression for mechanical behavior from continuous to discontinuous and from micro to macro, has the function of solving discontinuous problems and has high applicability and reliability in analyzing discontinuous and multi-scale problems. It can be used to study the deformation, damage, fracture and instability of homogeneous and heterogeneous target bodies, and simulate the whole process from crack initiation, expansion to structure failure. At present, near-field dynamics is rarely used in the field of land subsidence, and a complete framework has not been formed. The existing technology has two shortcomings: the numerical solution strategy of surface deformation is explicit time integration, which is difficult to provide stable solutions for land subsidence-ground fissure evolution problems which are regarded as quasi-static conditions; and the Mohr-Coulomb failure criterion based on strain needs to know a variety of geomechanical parameters of the case area in advance, and does not have universality. SUMMARY

[0006] In view of the above technical problems, the technical scheme adopted by the present application is:

[0007] According to a first aspect of the present application, a numerical simulation method for surface deformation evolution based on near-field dynamics is provided, and the method comprises the following steps:

[0008] In S100, a geological body model to be calculated is discretized into n1 material points, corresponding material point information is obtained, a virtual boundary layer with a preset thickness h is arranged outside the geological body model, the virtual boundary layer is discretized into n2 material points, and an initial geological body model containing n1+n2 material points is obtained; the material point information at least includes the number, coordinates and material parameters of the material points.

[0009] In S200, corresponding boundary conditions are applied to the virtual boundary layer, and the neighborhood of each material point in the initial geological body model is obtained, wherein the distance between any material point and any neighborhood material point in the corresponding neighborhood of the material point is less than a set distance d0.

[0010] In S300, the internal force acting on the material point at the current time step is obtained based on the change of the pore water pressure of the material point at the current time step, as the force acting on the material point at the current time step.

[0011] In S400, the displacement of the material point at the current time step is obtained by using an adaptive dynamic relaxation method based on the force acting on the material point at the current time step.

[0012] In S500, the stress of the material point at the current time step is obtained based on the displacement of the material point at the current time step.

[0013] In S600, the fracture of the virtual key is judged based on the obtained stress of the material point at the current time step, and the local damage value of the material point at the current time step is obtained; if the current time step is equal to a preset time step threshold, the current control program is exited, otherwise, S300 is executed.

[0014] According to a second aspect of the present application, a numerical simulation device for surface deformation evolution based on near-field dynamics is provided, and the device comprises:

[0015] An initial geological body model acquisition module is configured to discretize a geological body model to be calculated into n1 material points, obtain corresponding material point information, arrange a virtual boundary layer with a preset thickness h outside the geological body model, discretize the virtual boundary layer into n2 material points, and obtain an initial geological body model containing n1+n2 material points; the material point information at least includes the number, coordinates and material parameters of the material points.

[0016] A material point neighborhood acquisition module is configured to apply corresponding boundary conditions to the virtual boundary layer, and obtain the neighborhood of each material point in the initial geological body model, wherein the distance between any material point and any neighborhood material point in the corresponding neighborhood of the material point is less than a set distance d0.

[0017] The stress obtaining module obtains the internal force of the material point at the current time step based on the change of the pore water pressure of the material point at the current time step, and takes the internal force as the stress of the material point at the current time step.

[0018] The displacement obtaining module obtains the displacement of the material point at the current time step based on the stress of the material point at the current time step by using the adaptive dynamic relaxation method.

[0019] The stress obtaining module obtains the stress of the material point at the current time step based on the displacement of the material point at the current time step.

[0020] The local damage value obtaining module judges the fracture of the virtual key by using the Mohr-Coulomb failure criterion based on the stress of the material point at the current time step, and obtains the local damage value of the material point at the current time step.

[0021] The present application has at least the following beneficial effects:

[0022] The method and device for simulating the evolution of ground surface deformation based on near-field dynamics provided by the present application can effectively simulate the evolution process of ground surface deformation caused by groundwater exploitation by discretizing the geological model into a series of material points with material physical and mechanical information in space, taking the force of groundwater on the soil as the power source of the model, establishing the integral form of the balance equation based on the non-local action idea, and combining the adaptive dynamic relaxation method and the Mohr-Coulomb failure criterion based on stress.

[0023] It should be understood that the content described in this part is not intended to identify the key or important features of the embodiments of the present application, nor is it used to limit the scope of the present application. Other features of the present application will become apparent from the following description. BRIEF DESCRIPTION OF DRAWINGS

[0024] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings needed in the embodiment description will be briefly introduced below. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can also be obtained from these drawings without creative labor.

[0025] Figure 1A flow chart of the numerical simulation method of surface deformation evolution based on near-field dynamics provided by the embodiment of the present application is shown in the figure.

[0026] Figure 2 A schematic diagram showing the neighborhood of a material point is shown in the figure.

[0027] Figure 3 A schematic diagram of the geological body model used in the first embodiment is shown in the figure.

[0028] Figure 4 A schematic diagram of the calculation results of the first embodiment is shown in the figure.

[0029] Figure 5 A schematic diagram of the comparison between the first embodiment and the finite element-interface element model is shown in the figure. Figure 6

[0030] A schematic diagram of the geological body model used in the second embodiment is shown in the figure. Figure 7

[0031] A schematic diagram of the calculation results of the second embodiment is shown in the figure. Figure 8

[0032] A schematic diagram of the comparison between the second embodiment and the finite element-interface element model is shown in the figure. Figure 9 Figure 10 A schematic diagram of the geological body model used in the third embodiment is shown in the figure.

[0033] Figure 11 A schematic diagram of the calculation results of the third embodiment is shown in the figure.

[0034] Figure 12 A schematic diagram of the comparison between the third embodiment and the finite element-interface element model is shown in the figure.

[0035] DETAILED DESCRIPTION Figure 13 Figure 14 The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the accompanying drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person skilled in the art without creative work fall within the scope of protection of the present application. Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. The terminology used in the description herein is for describing the specific embodiments only and is not intended to be limiting of the application. As used herein, the term "and / or" includes any and all combinations of one or more of the associated listed items.

[0036]

[0037] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. The terminology used in the description herein is for describing the specific embodiments only and is not intended to be limiting of the application. As used herein, the term "and / or" includes any and all combinations of one or more of the associated listed items.

[0038] ​​It is to be understood that some of the example embodiments are described in terms of a process or method depicted as a flowchart. Although a flowchart can describe operations as a sequential process, many of the operations can be performed in parallel, concurrently or simultaneously. In addition, the order of the operations can be re-arranged. A process is terminated when its operations are completed, but could also terminate without waiting for the completion of all operations. A process can correspond to a method, a function, a procedure, a subroutine, a subprogram, etc.

[0039] The embodiment of the present application provides a surface deformation evolution numerical simulation method based on near-field dynamics, which comprises the following steps. Figure 1 As shown in the figure, the method can comprise the following steps:

[0040] S100, discretize a geological body model to be calculated into n1 material points, obtain corresponding material point information, set a virtual boundary layer with a preset thickness h outside the geological body model, discretize the virtual boundary layer into n2 material points, and obtain an initial geological body model containing n1+n2 material points; the material point information at least includes the number, coordinates and material parameters of the material points.

[0041] In the embodiment of the present application, the geological body model to be calculated can be discretized into a series of material points with material physics information in space. The material point refers to a node representing the volume and mass of soil material in a certain space range. The node itself does not occupy a space volume, and only the volume and mass of the material in the space range represented by the node are recorded in the material parameter matrix, so each material point can be regarded as the average value of the material physics and mechanics properties in a certain space range.

[0042] In the embodiment of the present application, the material parameters can include density, Young's modulus, Poisson's ratio, cohesion, internal friction angle, void ratio, bulk modulus, shear modulus, etc.

[0043] In the embodiment of the present application, the virtual boundary layer can include a left virtual boundary layer and a right virtual boundary layer located on the left and right sides of the geological body model, a front virtual boundary layer and a rear virtual boundary layer located on the front and rear sides of the geological body model, and a lower virtual boundary layer located on the lower side of the geological body model.

[0044] In the embodiment of the present application, the spacing △d between two adjacent material points in the initial geological body model can be set based on actual needs. In one illustrative embodiment, △d=5m.

[0045] In the embodiment of the present application, the material point information of each material point can be stored using a matrix. The virtual boundary layer has the same discretization manner as the geological body model to be calculated.

[0046] S200, applying corresponding boundary conditions on the virtual boundary layer and obtaining a neighborhood of each material point in the initial geological body model, wherein the distance between any material point and any neighborhood material point in the neighborhood of the material point is less than a set distance d0.

[0047] In the embodiment of the present application, the geostress of the geological body model to be calculated can be equivalent to the stress boundary condition of the calculation model, and the displacement constraint can be converted into the displacement boundary condition and applied on the virtual boundary layer and stored in a matrix. Specifically, the corresponding boundary conditions applied on the virtual boundary layer include horizontal fixed constraints applied on the left and right virtual boundary layers, longitudinal fixed constraints applied on the front and rear virtual boundary layers, and horizontal fixed constraints, longitudinal fixed constraints and vertical fixed constraints applied on the lower virtual boundary layer. The horizontal fixed constraint is used to constrain the model from moving along the X-axis, the longitudinal fixed constraint is used to constrain the model from moving along the Y-axis, and the vertical fixed constraint is used to constrain the model from moving along the Z-axis. In the embodiment of the present application, the coordinate system corresponding to the XYZ axes can be the coordinate system shown in FIG. 1. Figures 3 to 14

[0048] In the embodiment of the present application, the neighborhood refers to the near-field range in which a certain material point interacts with other material points, that is, the maximum distance d0 between the material point and other material points with which the material point interacts is a circle (two-dimensional) or a sphere (three-dimensional) with the radius d0, as shown in FIG. 2. Figure 2

[0049] In the embodiment of the present application, the interaction relationship between any material point and any neighborhood material point in the neighborhood of the material point is represented by a virtual bond.

[0050] In the embodiment of the present application, the thickness h of the virtual boundary layer can be equal to the radius of the neighborhood, that is, can be equal to d0.

[0051] S300, obtaining the internal force suffered by the material point at the current time step based on the change of the pore water pressure of the material point at the current time step, as the force suffered by the material point at the current time step.

[0052] In the embodiment of the present application, the initial model of the current geological body model is the initial geological body model.

[0053] In the embodiment of the present application, the internal force F i suffered by any material point i satisfies the following condition: F i =∫ h(i) j＝1 (F ij -F​​ji ) x dv j .

[0054] where F ij is the force of the jth neighbor particle of the neighborhood corresponding to the particle i on the particle i, F ji is the force of the particle i on the jth neighbor particle of the neighborhood of the particle i, i is valued from 1 to (n1+n2), j is valued from 1 to h(i), h(i) is the number of neighbor particles in the neighborhood corresponding to the particle i, v j is the volume of the jth neighbor particle, dv j represents the integration on v j .

[0055] where F ji satisfies the following conditions:

[0056] F ji = ((-3) x p i / wv i ) x w ji x ||x j0 -x i0 || + (15 x SM i ) / v i x w ji x e ij d ;

[0057] where w ji is the influence weight of the jth neighbor particle on the particle i. In an illustrative embodiment, w ji may be determined based on the distance between the jth neighbor particle and the particle i, the closer the distance, the greater the w ij . In another illustrative embodiment, the influence weight of each neighbor particle on the center point can be the same, for example, all 1. It should be noted that if the virtual bond between the neighbor particle and the center point is in a broken state, the corresponding influence weight is 0.

[0058] where |||| represents the norm, x i0 is the initial coordinate of the particle i, i.e., the coordinate in the initial geological model, and x j0 is the initial coordinate of the jth neighbor particle.

[0059] SM i is the shear modulus of the particle i, e ij d is the deviation component of the elongation state e ij of the virtual bond between the particle i and the jth neighbor particle. wv i is the weighted volume of the particle i, p iLet p be the near-field dynamic pressure of matter point i. i =-k i ×θ i +γ i ×pf i k i Let γ be the bulk modulus of point i. i Let pf be the fluid pressure coefficient at substance point i. i Let θ be the change in pore water pressure at material point i. i Let be the volume expansion rate of substance point i.

[0060] In this embodiment of the invention, the shear modulus and bulk modulus of the material points can be obtained based on the Young's modulus and Poisson's ratio of the material points, and the specific acquisition method can adopt existing methods. The fluid pressure coefficient of the material points can be a preset fluid pressure coefficient. In an illustrative embodiment, the fluid pressure coefficient of all material points can be set to 1.

[0061] In this embodiment of the invention, the change in pore water pressure at a material point is relative to the change in pore water pressure in the initial state. The change in pore water pressure in the initial state is 0. The change in pore pressure at each time step of the material point can be calculated by dividing the total change in pore water pressure by the total number of time steps and then multiplying by the current time step, i.e., Δpf=(Δpfs) / Ts×tc, where Δpf is the current change in pore pressure at the material point, Δpfs is the total change in pore water pressure, Ts is the total number of time steps (i.e., the preset time step threshold), and tc is the current time step of the material point.

[0062] Furthermore, in this embodiment of the invention, e ij =||(x j0 +x jc )-(x i0 +x ic )||-||x j0 -x i0 ||,e ij d =e ij -(θ i ×||x j0 -x i0 ||) / 3.

[0063] Furthermore, in this embodiment of the invention, wv i =∫1 h(i) w ij ×||x j0 -x i0 ||×||x j0 -x i0 ||×dv j ;

[0064] Further, in the embodiment of the present application, θ i = (3 / wv i ) x ∫1 h(i) w ij x ||x j0 -x i0 || x e ij x dv j .

[0065] wherein F ij satisfies the following condition:

[0066] F ij = ((-3) x p j / wv j ) x w ij x ||x j0 -x i0 || + (15 x SM j ) / wv j x w ij x e ij d ;

[0067] wherein w ij is the influence weight of the jth neighborhood material point on the material point i, SM j is the shear modulus of the jth neighborhood material point, wv j is the weighted volume of the jth neighborhood material point, p j is the near-field dynamic pressure of the jth neighborhood material point, p j = -k j x θ j + γ j x pf j , k j is the bulk modulus of the jth neighborhood material point, γ j is the fluid pressure coefficient of the jth neighborhood material point, pf j is the pore water pressure change value of the jth neighborhood material point, and θ j is the volumetric expansion rate of the jth neighborhood material point.

[0068] Further, in the embodiment of the present application, wv j = ∫1 g(j) w js x ||x j0 -x i0 || x ||x j0 -x i0 || x dv s . w js is the influence weight of the s th neighborhood material point in the neighborhood corresponding to the jth neighborhood material point on the jth neighborhood material point, and v sdv s is the volume of the s-th neighbor particle s is integrated, g(j) is the number of neighbor particles corresponding to the j-th neighbor particle.

[0069] Further, in the embodiment of the present application, θ j = (3 / wv j ) x ∫1 g(j) w js x ||x j0 -x i0 || x e ij x dv s .

[0070] S400, based on the force of the material point at the current time step, the displacement of the material point at the current time step is obtained by using the adaptive dynamic relaxation method.

[0071] In the embodiment of the present application, the adaptive dynamic relaxation method can be used to convert the near-field dynamic equilibrium equation into a balance equation in the form of ordinary differential equation by setting virtual damping and virtual mass, and the displacement of the material point is solved.

[0072] Further, in S400, the displacement x u of the material point at the u-th time step satisfies the following conditions:

[0073] x u = x u-1 + v u+(1 / 2) x △t;

[0074] wherein x u-1 is the displacement at the u-1-th time step, △t is the time step length, i.e. the step length between adjacent two time steps, v u+(1 / 2) is the velocity at the u+(1 / 2)-th time step, v u+(1 / 2) = ((2-dc u x △t) x v u-(1 / 2) + 2 x △t x DR -1 x F u ) / (2+dc u x △t), dc u is the virtual damping matrix at the u-th time step, F u is the internal force of the material point at the u-th time step, DR is the virtual diagonal density matrix, and v u-(1 / 2) is the velocity at the u-(1 / 2)-th time step. In an illustrative embodiment, △t can be set to 1.

[0075] S500, based on the displacement of the material point at the current time step, the stress of the material point at the current time step is obtained.

[0076] In S500, the stress tensor Si satisfy the following conditions:

[0077] S i = D: E i ;

[0078] wherein E i is the strain tensor of the material point i, E i = (1 / 2) (H i T H i -I), H i is a deformation gradient tensor obtained based on the displacement of the material point i at the current time step, H i = R·K -1 , K is a shape tensor before deformation of the material point i, denotes a tensor product, R is a shape tensor after deformation of the material point i, I is a unit matrix, D is a stiffness matrix, : denotes a double dot product, and · denotes a dot product.

[0079] In the embodiment of the present application, the stress tensor can include principal stresses acting on the x-axis, the y-axis and the z-axis and shear stresses acting on the xy, yz and zx planes, so that the stress of the material point including the principal stresses and the shear stresses can be obtained through the calculated stress tensor of the material point.

[0080] S600, based on the obtained stress of the material point at the current time step, judging the fracture of the virtual key by using the Mohr-Coulomb failure criterion to obtain the local damage value of the material point at the current time step; if the current time step is equal to the preset time step threshold, exiting the current control program, otherwise, executing S300.

[0081] Further, in the embodiment of the present application, the fracture of the virtual key can be judged based on the stress of the material point at the current time step, and the local damage value of the material point at the current time step can be obtained by using the Mohr-Coulomb failure criterion, which can specifically include the following steps:

[0082] S601, obtaining the stress tensor σ j suffered by the virtual key K b ij = (S i +S ij ) / 2, S i is the stress tensor of the material point i, and S ij is the stress tensor of the jth neighbor material point of the material point i.

[0083] S602, obtaining the normal stress σ j of K Nij =σ b ij ·n, and obtaining K j shear stress U ij =σ b ij -σ N ij ·n, where n is the direction vector.

[0084] S603, if ||U ij ||≥τ ij Or σ N ij ≥0, set K j The corresponding bond constant g ij =1, otherwise, that is, if ||U ij ||<τ ij And σ N ij <0, set g ij =0.

[0085] Among them, U ij For K j The corresponding shear stress, τ ij For K j shear strength, τ ij =σ N ij ×tan(φ ij )+c ij φ ij For K j The average of the internal friction angles of the two corresponding material points, c ij For K j The average cohesive force between the two corresponding material points.

[0086] In this embodiment of the invention, the bond constant is used to represent the integrity of the bond. If the bond constant is 1, it indicates that the bond is in a broken state. At this time, there is no longer an interaction relationship between the two interacting material points. That is, in this embodiment of the invention, when the shear stress of the bond exceeds the shear strength or the normal stress exceeds 0, the bond breaks. If the bond constant is 0, it indicates that the bond is in an unbroken state.

[0087] S604, Obtain the local damage value P of material point i at the current time step. i =1-(∫1 h(i) g ij ×dv j / ∫1 h(i) dv j ), dv j Let v be the volume of the j-th neighboring material point. j Integrate the points.

[0088] S605, if P i =0, indicating that material point i is in an unbroken state, meaning all virtual bonds are intact. If P i =1, indicating that matter point i is in a completely broken state, meaning all virtual bonds are broken, if 0 < P i <1 indicates that material point i is in a state of local destruction, i.e., some virtual bonds are destroyed.

[0089] In this embodiment of the invention, the preset time step threshold can be related to various factors such as the distance between material points and the size of the geological body model. Generally, it is necessary to test the impact of different time step thresholds on the stability of the results in advance. For example, the time step threshold can be set to 50, 200, 500, 1000, and 2000 respectively, and then the maximum displacement obtained by simulation can be compared. If the maximum displacement keeps changing when the threshold is between 50 and 500, while the displacements obtained at 1000 and 2000 are equal, then it proves that the time step threshold can be set to 1000. This way, it will not consume more model running time and can ensure the stability of the model. In an illustrative embodiment, the time step threshold can be set to 1000.

[0090] In this embodiment of the invention, the displacement of the material point corresponds to the soil settlement, the velocity of the material point corresponds to the soil settlement velocity, and the local damage of the material point corresponds to the local damage of the soil, i.e., ground fissures. The numerical simulation method for surface deformation evolution based on near-field dynamics provided by this invention can reveal the physical mechanism by which ground subsidence disasters transform into ground fissures.

[0091] (Example 1)

[0092] This embodiment provides a simulation of surface deformation caused by groundwater extraction under bedrock structural conditions, including the following steps:

[0093] (1) Model discretization and parameter initialization:

[0094] In this embodiment, as Figure 3 As shown, the geological model has a length of 2000m, a height of 500m, a thickness of 100m, and a bedrock burial depth of 100m. It is divided into 400, 100, and 20 material points along the length, height, and thickness directions, respectively. The virtual boundary layer is divided into 3 material points, with a spacing of 5m between the material points. The near-field range, or neighborhood, is 3.015 times the spacing between the material points.

[0095] Initialize the material parameter matrices as shown in Table 1:

[0096] Table 1:

[0097]

[0098] (2) Neighborhood point retrieval:

[0099] According to the given near-field range, the number of other material points in the space range of each material point is retrieved and stored in the neighborhood matrix.

[0100] (3) Boundary condition application:

[0101] The left and right virtual boundary layers of the model are applied with horizontal fixed constraints, the front and back virtual boundary layers are applied with longitudinal fixed constraints, and the lower virtual boundary layer is applied with horizontal, longitudinal, and vertical fixed constraints. In this embodiment, since the rock is very hard and generally does not displace due to groundwater extraction, the rock area can be treated as a boundary, i.e., direct displacement constraints are applied, i.e., horizontal, longitudinal, and vertical fixed constraints are applied to the rock area.

[0102] (4) Time domain integration:

[0103] In this embodiment, the groundwater level change of the aquifer is 100 m. The groundwater level change is used as the power source of the model to solve the force of the material points at each time step; the adaptive dynamic relaxation algorithm is used to dynamically solve the virtual mass density matrix and the virtual damping coefficient to solve the displacement of the material points at each time step; based on the displacement, the stress of the material points at each time step is solved, and then the stress of each virtual key is calculated. In this embodiment, the time step threshold is set to 1000 steps.

[0104] (5) Damage judgment:

[0105] The Mohr-Coulomb failure criterion based on stress is used to judge the fracture of the virtual key of the material point. If the shear stress exceeds the shear strength or the normal stress is greater than 0, the key constant of the corresponding key is 0, otherwise it is 1, and the local damage value of each material point is obtained by integration.

[0106] (6) Result analysis:

[0107] After the calculation is completed, the evolution process of the differential settlement disaster caused by groundwater exploitation under the bedrock structure condition is obtained, as shown in FIG. 4. Figure 4 It can be seen that under the action of groundwater exploitation, the soil is compressed, the upper soil of the bedrock gradually breaks, and finally a ground fissure is formed.

[0108] To verify the effectiveness of the proposed near-field dynamics-based surface deformation numerical simulation method and device, a finite element-interface element model is constructed for comparison. The maximum and minimum values of horizontal and vertical deformation are compared as shown in Table 2, and the spatial distribution is compared as shown in FIGS. 5 and 6. Figure 5 Figure 6

[0109] Table 2:

[0110] ​​

[0111] Example 2:

[0112] This embodiment provides a simulation of surface deformation caused by groundwater extraction under fractured tectonic conditions, including the following steps:

[0113] (1) Model discretization and parameter initialization:

[0114] In this embodiment, as Figure 7 As shown, the geological model is 2000m long, 500m high, and 100m thick. The fault structure is a normal fault with a dip angle of 100°. It is divided into 400, 100, and 20 material points along the length, height, and thickness directions, respectively. The virtual boundary layer is divided into 3 material points. The distance between the material points is 5m, and the near-field range is 3.015 times the distance between the material points.

[0115] Initialize each material parameter matrix. The material parameter matrices in this embodiment are the same as those in Embodiment 1.

[0116] (2) Neighborhood point retrieval:

[0117] Based on the given near-field range, retrieve the numbers of other matter points within that spatial range for each matter point and store them in the neighborhood matrix.

[0118] (3) Apply boundary conditions:

[0119] Apply horizontal fixed constraints to the left and right virtual boundary layers of the model, apply vertical fixed constraints to the front and rear virtual boundary layers, and apply horizontal, vertical, and vertical fixed constraints to the lower virtual boundary layer.

[0120] (4) Time-domain integral:

[0121] In this implementation, the groundwater level change in the aquifer is 100m. The groundwater level change is used as the dynamic source of the model to solve for the force on the material point at each time step. An adaptive dynamic relaxation algorithm is used to dynamically solve for the virtual mass density matrix and virtual damping coefficient, and to solve for the displacement of the material point at each time step. Based on the displacement, the stress on the material point at each time step is solved, and then the stress on each virtual bond is calculated. In this embodiment, the time step threshold is set to 1000 steps.

[0122] (5) Damage assessment:

[0123] The fracture status of virtual bonds at material points is determined using a stress-based Mohr-Coulomb failure criterion. If the shear stress exceeds the shear strength or the normal stress is greater than 0, the bond constant of the corresponding bond is 0; otherwise, it is 1. The local damage value of each material point is obtained by integration.

[0124] (6) Result analysis:

[0125] After the calculation, the evolution process of the differential settlement disaster caused by groundwater exploitation under the bedrock structure condition is obtained, as shown in Figure 8 It can be seen that under the action of groundwater exploitation, the soil in the lower plate is compressed, and the soil in the upper plate near the fracture zone gradually breaks, and finally forms a ground fissure.

[0126] To verify the effectiveness of the proposed surface deformation numerical simulation method and device based on near-field dynamics, a finite element-interface element model is constructed for comparison. The maximum values of horizontal and vertical deformation are shown in Table 3, and the spatial distribution is shown in Figure 9 and Figure 10 .

[0127] Table 3:

[0128]

[0129]

[0130] Example Three:

[0131] The present embodiment provides a simulation of surface deformation caused by groundwater exploitation under outcrop structure conditions, comprising the following steps:

[0132] (1) Model discretization:

[0133] In this embodiment, as shown in Figure 11 , the length of the geological body model is 2000m, the height is 500m, the thickness is 100m, and the outcrop dip angles are 75° and 105° respectively. Along the length, height and thickness directions, 400, 100 and 20 material points are divided respectively, the virtual boundary layer is divided into 3 material points, the material point spacing is 5m, and the near-field range is 3.015 times the material point spacing.

[0134] Initialize the material parameter matrix, and the material parameter matrix of the present embodiment is the same as that of example one.

[0135] (2) Neighborhood point retrieval:

[0136] According to the given near-field range, the other material point numbers within the space range of each material point are retrieved and stored in the neighborhood matrix.

[0137] (3) Boundary condition application:

[0138] The left and right virtual boundary layers of the model are applied with horizontal fixed constraints, the front and rear virtual boundary layers are applied with longitudinal fixed constraints, the lower virtual boundary layer is applied with horizontal fixed constraints, longitudinal fixed constraints and vertical fixed constraints, and the rock area is applied with horizontal fixed constraints, longitudinal fixed constraints and vertical fixed constraints.

[0139] (4) Time domain integration:

[0140] In this embodiment, the water table of the aquifer changes by 100 m. The change of the water table is taken as the driving source of the model to solve the force of the material point at each time step; the self-adaptive dynamic relaxation algorithm is used to dynamically solve the virtual mass density matrix and the virtual damping coefficient to solve the displacement of the material point at each time step; based on the displacement, the stress of the material point at each time step is solved, and then the stress borne by each virtual key is calculated. In this embodiment, the time step threshold is set to 1000 steps.

[0141] (5) Damage judgment:

[0142] The fracture of the virtual key of the material point is judged by the Mohr-Coulomb failure criterion based on stress. If the shear stress exceeds the shear strength or the normal stress is greater than 0, the key constant of the corresponding key is 0, otherwise it is 1, and the local damage value of each material point is obtained by integration.

[0143] (6) Result analysis:

[0144] After the calculation is completed, the evolution process of the differential settlement disaster caused by groundwater exploitation under the bedrock structure condition is obtained, as shown in FIG. 6. Figure 12 It can be seen that under the action of groundwater exploitation, the soil is compressed, the soil near the bedrock gradually breaks, and finally the ground fissure is formed.

[0145] To verify the effectiveness of the numerical simulation method and system of surface deformation based on near-field dynamics proposed, a finite element-interface element model is constructed for comparison. The maximum and minimum values of horizontal and vertical deformation are shown in Table 4, and the spatial distribution is shown in FIGS. 7 and 8. Figure 13 Figure 14

[0146] Table 4:

[0147]

[0148] The above three embodiments all achieve simulation effects close to the traditional finite element model. It can be seen that the numerical simulation method of surface deformation evolution based on near-field dynamics can effectively simulate the process of differential ground settlement disaster caused by groundwater exploitation under various geological structure conditions such as bedrock, fracture and outcrop, and at the same time, avoids the defect that the finite element-interface element model needs to know in advance whether the crack exists and its position and size.

[0149] Based on the same inventive concept, the embodiment of the present application provides a numerical simulation device of surface deformation evolution based on near-field dynamics, which comprises:

[0150] ​​An initial geologic body model acquisition module is configured to discretize a geologic body model to be calculated into n1 material points, to obtain corresponding material point information, and to set a virtual boundary layer with a preset thickness h outside the geologic body model and discretize the virtual boundary layer into n2 material points to obtain an initial geologic body model containing n1+n2 material points; the material point information at least includes the number, coordinates and material parameters of the material points.

[0151] A material point neighborhood acquisition module is configured to apply corresponding boundary conditions on the virtual boundary layer and to obtain a neighborhood of each material point in the initial geologic body model, wherein the distance between any material point and any neighborhood material point in the neighborhood of the material point is less than a set distance d0.

[0152] A material point force acquisition module is configured to obtain, as the force on a material point at a current time step, the internal force on the material point at the current time step based on the change in pore water pressure of the material point at the current time step.

[0153] A material point displacement acquisition module is configured to obtain, by using an adaptive dynamic relaxation method, the displacement of a material point at a current time step based on the force on the material point at the current time step.

[0154] A material point stress acquisition module is configured to obtain the stress of a material point at a current time step based on the displacement of the material point at the current time step.

[0155] A material point damage judgment module is configured to obtain, by using the Mohr-Coulomb failure criterion, the local damage value of a material point at a current time step based on the stress of the material point at the current time step, and to judge the fracture of a virtual key. Figure 1 The functions and the like that can be achieved by the functional modules of the device can refer to the descriptions of the embodiments shown in Figure 1 The descriptions of the embodiments shown in

[0156] The embodiments of the present application also provide an electronic device, which includes at least one processor and a memory in communication connection with the at least one processor; the memory stores instructions executable by the at least one processor, and the instructions are configured to execute the method described in the embodiments of the present application.

[0157] The embodiments of the present application also provide a non-transitory computer readable storage medium storing computer executable instructions, and the computer executable instructions are configured to execute the method described in the embodiments of the present application.

[0158] It should be understood that the various forms of flow shown above can be re-ordered, added to, or have steps deleted, for example. The steps recited in the present disclosure can be performed in parallel, in series, or in a different order, as long as the desired results of the technology disclosed herein are achieved, and this document is not limited in this regard.

[0159] The specific embodiments discussed hereinabove are illustrative of various aspects of the present application. Alterations, modifications, combinations, sub-combinations, and the like are intended to be within the scope of the present application. Accordingly, although specific embodiments have been described, these are examples only and are not limiting upon the scope of the present application.

Claims

1. A numerical simulation method of surface deformation evolution based on near-field dynamics, characterized in that, The method comprises the following steps: S100, discretizing a geological body model to be calculated into n1 material points to obtain corresponding material point information, and setting a virtual boundary layer with a preset thickness h outside the geological body model and discretizing the virtual boundary layer into n2 material points to obtain an initial geological body model containing n1+n2 material points; the material point information at least comprises a number, coordinates and material parameters of the material points; S200, applying corresponding boundary conditions on the virtual boundary layer and obtaining a neighborhood of each material point in the initial geological body model, wherein the distance between any material point and any neighborhood material point in the neighborhood of the material point is less than a set distance d0; S300, obtaining an internal force suffered by the material point at a current time step based on the change of the pore water pressure of the material point at the current time step as the force suffered by the material point at the current time step; S400, obtaining the displacement of the material point at the current time step by using an adaptive dynamic relaxation method based on the force suffered by the material point at the current time step; S500, obtaining the stress of the material point at the current time step based on the displacement of the material point at the current time step; S600, judging the fracture of a virtual key by using a Mohr-Coulomb failure criterion based on the stress of the material point at the current time step to obtain a local damage value of the material point at the current time step; if the current time step is equal to a preset time step threshold, exiting the current control program, otherwise, executing S300; In S400, the displacement x of the material point at the u-th time step u satisfies the following condition: x u =x u-1 +v u+（1 / 2） ×△t; wherein x u-1 is the displacement at the u-1 time step, At is the time step, v u+（1 / 2） is the velocity at the u+(1 / 2) time step, v u +（1 / 2） = ((2-dc u x At) x v u-（1 / 2） + 2 x At x DR -1 x F u ) / (2 + dc u x At), dc u is the virtual damping matrix at the u time step, F u is the internal force experienced by the material point at the u time step, DR is the virtual diagonal density matrix, v u-（1 / 2） is the velocity at the u-(1 / 2) time step; The step of judging the fracture of the virtual key by using the Mohr-Coulomb failure criterion based on the stress of the material point at the current time step to obtain the local damage value of the material point at the current time step comprises the following steps: S601, obtaining a virtual key K corresponding to a material point i j The stress tensor σ received b ij = (S i + S ij ) / 2, S i is the stress tensor of the material point i, S ij is the stress tensor of the j-th neighborhood material point of the material point i, and j takes values from 1 to h(i); S602, obtaining K j Normal stress σ N ij =σ b ij •n, and obtaining K j Shear stress U ij =σ b ij -σ N ij •n, n being a direction vector; S603, if ||U ij ||≥τ ij or σ N ij ≥0, set K j =1, if ||U ij ||<τ ij and σ ij N <0, set g ij =0, U ij is the shear stress corresponding to the virtual key j, τ ij is the shear strength of K ij , τ j =σ ij N ij ×tan(φ ij )+c ij , φ ij is the average of the internal friction angle of the two material points corresponding to K j , c ij is the average of the cohesion of the two material points corresponding to K j ;​ S604, obtaining the local damage value P of the material point i i =1-(∫1 h（i） g ij ×dv j / ∫1 h（i） dv j ); S605, if P i = 0, determine that the material point i is in an undamaged state, if P i = 1, determine that the material point i is in a completely damaged state, if 0 < P i < 1, determine that the material point i is in a partially damaged state.

2. The method of claim 1, wherein, d0=f×△d, f is a preset value greater than 1, and △d is the spacing between two adjacent material points in the initial geological body model.

3. The method of claim 1, wherein, The interaction relationship between any material point and any neighborhood material point in the neighborhood of the material point is represented by a virtual key; The internal force F received by any material point i i Satisfies the following condition: F i =∫ h（i） j=1 (F ij -F ji )×dv j ; wherein F ij is the force of the jth neighbor material point within the neighborhood corresponding to material point i on material point i, F ji is the force of the jth neighbor material point within the neighborhood of material point i, i has a value from 1 to (n1+n2), j has a value from 1 to h(i), h(i) is the number of neighbor material points within the neighborhood corresponding to material point i, v j is the volume of the jth neighbor material point, dv j denotes integration over v j ; F ji satisfies the following conditions: F ji = ((-3) x p i / wv i ) x w ji x ||x j0 - x i0 || + (15 x SM i ) / wv i x w ji x e ij d ; w ji is the influence weight of the jth neighborhood material point on the material point i, || || represents the norm, x i0 is the initial coordinate of the material point i, that is, the coordinate in the initial geological body model, x j0 is the initial coordinate of the jth neighborhood material point, SM i is the shear modulus of the material point i, e ij d is the elongation state of the virtual bond between the material point i and the jth neighborhood material point e ij is the deviation component of wv i is the weighted volume of the material point i, p i is the near-field dynamics pressure of the material point i, p i =-k i ×θ i +γ i ×pf i , k i is the bulk modulus of the material point i, γ i is the fluid pressure coefficient of the material point i, pf i is the pore water pressure change value of the material point i, θ i is the volume expansion rate of the material point i; F ij satisfies the following conditions: F ij = ((-3) x p j / wv j ) x w ij x ||x j0 - x i0 || + (15 x SM j ) / wv j x w ij x e ij d ; where w ij is the influence weight of the jth neighbor material point on material point i, SM j is the shear modulus of the jth neighbor material point, wv j is the weighted volume of the jth neighbor material point, p j is the near-field dynamic pressure of the jth neighbor material point, p j = -k j x θ j + γ j x pf j , k j is the bulk modulus of the jth neighbor material point, γ j is the fluid pressure coefficient of the jth neighbor material point, pf j is the change in pore water pressure of the jth neighbor material point, θ j is the volumetric swelling ratio of the jth neighbor material point.

4. The method of claim 3, wherein, e ij =||(x j0 +x jc )-(x i0 +x ic )||-||x j0 -x i0 ||,e ij d =e ij -(θ i ×||x j0 -x i0 ||) / 3,where, x ic Let x be the displacement of point i. jc Let be the displacement of the j-th neighboring material point; wv i =∫1 h（i） w ij ×||x j0 -x i0 ||×||x j0 -x i0 ||×dv j ; θ i = (3 / wv i ) x ∫1 h（i） w ij x ||x j0 -x i0 || • e ij x dv j .

5. The method of claim 3, wherein, In S500, the stress tensor S of any material point i i satisfies the following condition: S i =D:E i ; E i is the strain tensor for material point i, E i = (1 / 2) (H i T H i -I), H i is the deformation gradient tensor obtained based on the displacement of material point i at the current time step, H i =R•K -1 , K is the shape tensor of material point i before deformation, K=∫1 h（i） w ij ×((x j0 -x i0 ) (x j0 -x i0 ))×dv j , denotes the tensor product, R is the shape tensor of material point i after deformation, R=∫1 h（i） w ij ×(((x j0 +x jc )- (x i0 +x ic ) (x j0 -x i0 ))×dv j , I is the unit matrix, D is the stiffness matrix, : denotes the double dot product, • denotes the dot product.

6. The method of claim 1, wherein, The virtual boundary layer comprises a left virtual boundary layer and a right virtual boundary layer located on the left and right sides of the geological body model, a front virtual boundary layer and a rear virtual boundary layer located on the front and rear sides of the geological body model, and a lower virtual boundary layer located on the lower side of the geological body model; The step of applying corresponding boundary conditions on the virtual boundary layer comprises applying horizontal fixed constraints on the left virtual boundary layer and the right virtual boundary layer, applying longitudinal fixed constraints on the front virtual boundary layer and the rear virtual boundary layer, and applying horizontal fixed constraints, longitudinal fixed constraints and vertical fixed constraints on the lower virtual boundary layer.

7. A device for numerical simulation of evolution of surface deformation based on near-field dynamics, characterized in that, The device comprises: An initial geological body model acquisition module, configured to discretize a geological body model to be calculated into n1 material points to obtain corresponding material point information, and set a virtual boundary layer with a preset thickness h outside the geological body model and discretize the virtual boundary layer into n2 material points to obtain an initial geological body model containing n1+n2 material points; the material point information at least comprises a number, coordinates and material parameters of the material points; An initial geological body model acquisition module, configured to discretize a geological body model to be calculated into n1 material points to obtain corresponding material point information, and set a virtual boundary layer with a preset thickness h outside the geological body model and discretize the virtual boundary layer into n2 material points to obtain an initial geological body model containing n1+n2 material points; the material point information at least comprises a number, coordinates and material parameters of the material points; a material point neighborhood acquisition module, configured to apply a corresponding boundary condition on the virtual boundary layer and acquire a neighborhood of each material point in the initial geologic body model, wherein a distance between any material point and any neighborhood material point in the neighborhood of the material point is less than a set distance d0; a material point force acquisition module, configured to acquire, as a force borne by the material point at a current time step, an internal force borne by the material point at the current time step based on a change in pore water pressure of the material point at the current time step; a material point displacement acquisition module, configured to acquire, by using an adaptive dynamic relaxation method, a displacement of the material point at the current time step based on the force borne by the material point at the current time step; a material point stress acquisition module, configured to acquire a stress of the material point at the current time step based on the displacement of the material point at the current time step; a material point damage judgment module, configured to acquire a local damage value of the material point at the current time step by using a Mohr-Coulomb failure criterion to judge a fracture condition of a virtual key based on the stress of the material point at the current time step; where the displacement of the material point at the u-th time step x u satisfies the following condition: x u =x u-1 +v u+（1 / 2） ×△t; wherein x u-1 is the displacement at the u-1 time step,△t is the time step length, v u+（1 / 2） is the velocity at the u+ (1 / 2) time step, v u +（1 / 2） = ((2-dc u x△t) x v u-（1 / 2） + 2 x△t x DR -1 x F u ) / (2 + dc u x△t), dc u is the virtual damping matrix at the u time step, F u is the internal force experienced by the material point at the u time step, DR is the virtual diagonal density matrix, v u-（1 / 2） is the velocity at the u- (1 / 2) time step; the material point damage judgment module is specifically configured to perform the following operations: S601, obtaining a virtual key K corresponding to a material point i j a stress tensor σ received b ij = (S i +S ij ) / 2, S i is a stress tensor of the material point i, S ij is a stress tensor of a j-th neighborhood material point of the material point i, and j is valued from 1 to h(i); S602, obtaining K j the normal stress σ N ij =σ b ij •n, and obtaining K j the shear stress U ij =σ b ij -σ N ij •n, n being a direction vector; S603, if ||U ij ||≥τ ij or σ N ij ≥0, set K j The corresponding key constant g ij =1, if ||U ij ||<τ ij and σ N ij <0, set g ij =0, U ij is the shear stress corresponding to the virtual key j, τ ij is the shear strength of K j , τ ij =σ N ij ×tan(φ ij ) + c ij , φ ij is the average of the internal friction angle of the two material points corresponding to K j , c ij is the average of the cohesion of the two material points corresponding to K j ; S604, obtaining the local damage value P of the material point i i =1-(∫1 h（i） g ij ×dv j / ∫1 h（i） dv j ); S605, if P i = 0, determine that the material point i is in an undamaged state, if P i = 1, determine that the material point i is in a completely damaged state, and if 0 < P i < 1, determine that the material point i is in a partially damaged state.

Images