Ground surface deformation evolution numerical simulation method and device based on near-field dynamics

Abstract

The invention provides an earth surface deformation evolution numerical simulation method based on near-field dynamics, which comprises the following steps: dispersing a geologic body model and a virtual boundary layer which need to be calculated into a plurality of material points to obtain an initial geologic body model; acquiring a neighborhood of each substance point in the initial geologic body model; based on the pore water pressure change of the substance point at the current time step, acquiring the internal force of the substance point at the current time step as the stress of the substance point at the current time step; based on the stress of the substance point at the current time step, obtaining the displacement of the substance point at the current time step by using an adaptive dynamic relaxation method; based on the displacement of the substance point in the current time step, obtaining the stress of the substance point in the current time step; and based on the stress of the substance point at the current time step, judging the fracture condition of the virtual key by utilizing a Mohr-Coulomb failure criterion, and obtaining a local damage value of the substance point at the current time step. The invention provides another scheme capable of revealing a physical mechanism that ground subsidence disasters into ground fissures.

Description

[0001] Technology Neighborhood

[0002] The invention relates to a geological disaster neighborhood, and in particular to a method and a device for numerical simulation of surface deformation evolution based on near-field dynamics. Background Art

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

[0004] Numerical simulation methods are an important means to study the evolution and development of ground subsidence and ground fissures. Existing ground fissure simulation methods are mainly based on the finite element method and finite difference method based on the theory of continuous medium mechanics. However, the process of ground subsidence disaster turning into ground fissures is a nonlinear discontinuous time-space problem. The above methods have a contradiction between continuous theory and discontinuous practice. In addition, when using the theory of discontinuous medium mechanics, such as the discrete element method, interface element method and extended finite element method, on the one hand, it is necessary to know in advance whether cracks exist and their location and size. On the other hand, after the cracks propagate, the grid must be re-divided or additional functions must be introduced to describe the discontinuous problem, making it difficult to achieve automatic and accurate simulation of complex cracks.

[0005] Peridynamics is based on the idea of ​​nonlocal action and describes the mechanical behavior of materials by solving spatial integral equations that can be applied to discontinuities. The peridynamic method has a unified description of mechanical behavior from continuous to discontinuous, microscopic to macroscopic, has meshless properties and the function of solving discontinuous problems, and shows 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 targets, and simulate the entire process of crack initiation, expansion and structural destruction. At present, the application of peridynamics in the neighborhood of ground subsidence is very rare, and a complete framework system has not yet been formed. In addition, the existing technology has two shortcomings: the numerical solution strategy for surface deformation is an explicit time integration, which is difficult to provide a stable solution for the ground subsidence-ground fissure evolution problem regarded as a quasi-static condition; the strain-based Mohr-Coulomb failure criterion requires a prior understanding of multiple geomechanical parameters of the case area and is not universal. Summary of the invention

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

[0007] According to a first aspect of the present invention, a method for numerical simulation of surface deformation evolution based on peridynamics is provided, the method comprising the following steps:

[0008] S100, discretize the 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 on the outside of 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 includes the number, coordinates and material parameters of the material point.

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

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

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

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

[0013] S600, judging the fracture condition of the virtual key based on the obtained stress of the material point at the current time step, and obtaining 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, exit the current control program, otherwise, execute S300.

[0014] According to a second aspect of the present invention, there is provided a device for numerically simulating the evolution of surface deformation based on peridynamics, the device comprising:

[0015] The initial geological body model acquisition module is used to discretize the 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 on the outside of 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 includes the number, coordinates and material parameters of the material point.

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

[0017] The material point force acquisition module is used to obtain the internal force received 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 received by the material point at the current time step.

[0018] The material point displacement acquisition module is used to obtain the displacement of the material point in the current time step based on the force on the material point in the current time step using the adaptive dynamic relaxation method.

[0019] The material point stress acquisition module is used to obtain 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 material point local damage value acquisition module is used to determine the fracture of the virtual bond based on the stress of the material point at the current time step using the Mohr-Coulomb failure criterion to obtain the local damage value of the material point at the current time step.

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

[0022] The embodiment of the present invention provides a method and device for numerical simulation of surface deformation evolution based on near-field dynamics. The method discretizes the geological body model into a series of material points in space with material physical and mechanical information, takes the force of groundwater on the soil as the power source of the model, establishes an integral form of equilibrium equation based on the idea of ​​non-local action, and combines the adaptive dynamic relaxation method and the stress-based Mohr-Coulomb failure criterion to achieve effective simulation of the evolution process of surface deformation caused by groundwater exploitation. The method can reveal the physical mechanism of the transformation of ground subsidence disasters into ground fissures, and at least solve the current problems of difficulty in automatically and accurately simulating the process of ground subsidence disasters into ground fissures based on continuous medium theory and discontinuous medium theory, and the problem of failure to reveal the evolution mechanism of ground subsidence-ground fissures, as well as the instability and non-universality existing in the current numerical simulation of ground subsidence-ground fissures based on near-field dynamics.

[0023] It should be understood that the contents described in this section are not intended to identify the key or important features of the embodiments of the present invention, nor are they intended to limit the scope of the present invention. Other features of the present invention will become easily understood through the following description. BRIEF DESCRIPTION OF THE DRAWINGS

[0024] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings required for use in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings in the following description are only some embodiments of the present invention. For ordinary technicians in this field, other accompanying drawings can be obtained based on these accompanying drawings without paying any creative work.

[0025] Figure 1A flow chart of a method for numerical simulation of surface deformation evolution based on peridynamics provided in an embodiment of the present invention;

[0026] Figure 2 A schematic diagram showing the neighborhood of a material point;

[0027] Figure 3 Schematic diagram of the geological body model used in Example 1;

[0028] Figure 4 It is a schematic diagram of calculation results of Example 1;

[0029] Figure 5 and Figure 6 It is a comparative schematic diagram of the first embodiment and the finite element-interface element model;

[0030] Figure 7 Schematic diagram of the geological body model used in Example 2;

[0031] Figure 8 It is a schematic diagram of the calculation results of Example 2;

[0032] Fig. 9 and Fig.10 A schematic diagram for comparing the second embodiment and the finite element-interface element model;

[0033] Fig.11 Schematic diagram of the geological body model used in Example 3;

[0034] Fig.12 It is a schematic diagram of the calculation results of Example 3;

[0035] Fig.13 and Fig.14 It is a schematic diagram for comparing the third embodiment and the finite element-interface element model. DETAILED DESCRIPTION

[0036] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by technicians in this field without creative work are within the scope of protection of the present invention.

[0037] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as those generally understood by those skilled in the art of the present invention. The terms used herein in the specification of the present invention are only for the purpose of describing specific embodiments and are not intended to limit the present invention. The term "and / or" used herein includes any and all combinations of one or more related listed items.

[0038] It should be noted that some exemplary embodiments are described as processes or methods depicted as flow charts. Although the flow charts describe the steps as sequential processes, many of the steps can be implemented in parallel, concurrently or simultaneously. In addition, the order of the steps can be rearranged. The process can be terminated when its operation is completed, but it can also have additional steps not included in the accompanying drawings. The process can correspond to a method, function, procedure, subroutine, subprogram, etc.

[0039] The embodiment of the present invention provides a method for numerical simulation of surface deformation evolution based on peridynamics, such as Figure 1 As shown, the method may include the following steps:

[0040] S100, discretize the 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 on the outside of 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 includes the number, coordinates and material parameters of the material point.

[0041] In the embodiment of the present invention, the geological body model to be calculated can be discretized into a series of material points in space with material physics information. Material points refer to nodes that represent the volume and mass of soil materials within a certain spatial range. The nodes themselves do not occupy the spatial volume, and only the material volume and mass within the spatial range they represent are recorded in the material parameter matrix. Therefore, each material point can be regarded as the mean value of the physical and mechanical properties of the material within a certain spatial range.

[0042] In an embodiment of the present invention, the material parameters may include density, Young's modulus, Poisson's ratio, cohesion, internal friction angle, void ratio, bulk modulus, shear modulus, and the like.

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

[0044] In the embodiment of the present invention, the distance Δd between two adjacent material points in the initial geological body model can be set based on actual needs. In an exemplary embodiment, Δd=5m.

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

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

[0047] In an embodiment of the present invention, the ground stress 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 a displacement boundary condition, applied to the virtual boundary layer, and stored in a matrix. Specifically, the corresponding boundary conditions applied to the virtual boundary layer include: applying horizontal fixed constraints to the left virtual boundary layer and the right virtual boundary layer, applying longitudinal fixed constraints to the front virtual boundary layer and the rear virtual boundary layer, and applying horizontal fixed constraints, longitudinal fixed constraints and vertical fixed constraints to the lower virtual boundary layer. Among them, the horizontal fixed constraint is used to constrain the model from moving in the X-axis, the longitudinal fixed constraint is used to constrain the model from moving in the Y-axis, and the vertical fixed constraint is used to constrain the model from moving in the Z-axis. In an embodiment of the present invention, the coordinate system corresponding to the XYZ axis can be Figures 3 to 14 The coordinate system shown.

[0048] In the embodiment of the present invention, the neighborhood refers to the near-field range of a material point interacting 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 a radius, such as Figure 2 As shown, the neighborhood of material point 1 and material point 2 is the area corresponding to the circle drawn with them as the center point and d0 as the radius. In the embodiment of the present invention, d0=f×△d, and f is a preset value greater than 1. In an exemplary embodiment, f=3.015.

[0049] In the embodiment of the present invention, the interaction relationship between any material point and any neighboring material point in its corresponding neighborhood is represented by a virtual key.

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

[0051] S300, based on the change of pore water pressure of the material point in the current time step, the internal force received by the material point in the current time step is obtained as the force received by the material point in the current time step.

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

[0053] In the embodiment of the present invention, the internal force F exerted on any material point i is i Meet the following conditions: F i =∫ h(i) j＝1 (F ij -Fji )×dv j .

[0054] Among them, F ij is the force exerted by the jth neighboring material point in the neighborhood corresponding to material point i on material point i, F ji is the force exerted by material point i on the jth neighboring material point in its neighborhood, i ranges from 1 to (n1+n2), j ranges from 1 to h(i), h(i) is the number of neighboring material points in the neighborhood corresponding to material point i, v j is the volume of the jth neighborhood material point, dv j Indicates v j Perform integration.

[0055] Among them, F ji The following conditions must be met:

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

[0057] Among them, w ji is the influence weight of material point i on the jth neighboring material point. ji It can be determined based on the distance between the jth neighboring material point and material point i. The closer the distance, the greater the value of w ij In another exemplary embodiment, the influence weight of each neighboring material point on its center point may be the same, for example, all 1. It should be noted that if the virtual bond between the neighboring material point and the center point is in a broken state, the corresponding influence weight is 0.

[0058] Among them, |||| represents the norm, x i0 is the initial coordinate of material point i, i.e., the coordinate in the initial geological model, x j0 is the initial coordinate of the jth neighborhood material point.

[0059] SM i is the shear modulus of material point i, e ij d is the elongation state e of the virtual bond between material point i and the jth neighboring material point ij The deviation component of wv i is the weighted volume of material point i, p iis the peridynamic pressure at material point i, p i =-k i ×θ i +γ i ×pf i , k i is the bulk modulus of material point i, γ i is the fluid pressure coefficient at material point i, pf i is the change in pore water pressure at material point i, θ i is the volume expansion rate of material point i.

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

[0061] In the embodiment of the present invention, the pore water pressure change value of the material point is the change value of the pore water pressure relative to the initial state. The change value of the pore water pressure in the initial state is 0, and the pore pressure change value of the material point at each time step can be calculated according to the total pore water pressure change value divided by the total time step and then multiplied by the current time step, that is, △pf=(△pfs) / Ts×tc, △pf is the current pore pressure change value of the material point, △pfs is the total pore water pressure change value, Ts is the total time step, that is, the preset time step threshold, and tc is the current time step of the material point.

[0062] Furthermore, in the embodiment of the present 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] Further, in the embodiment of the present invention, wv i =∫1 h(i) w ij ×||x j0 -x i0 ||×||x j0 -x i0 ||×dv j ;

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

[0065] Among them, F ij The following conditions must be met:

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

[0067] Among them, w ij is the influence weight of the jth neighboring material point on material point i, SM j is the shear modulus of the jth neighboring material point, wv j is the weighted volume of the jth neighborhood material point, p j is the peridynamic pressure of the jth neighboring material point, p j =-k j ×θ j +γ j ×pf j , k j is the bulk modulus of the jth neighborhood material point, γ j is the fluid pressure coefficient of the jth neighboring material point, pf j is the change in pore water pressure of the jth neighboring material point, θ j is the volume expansion rate of the jth neighborhood material point.

[0068] Further, in the embodiment of the present invention, wv j =∫1 g(j) w js ×||x j0 -x i0 ||×||x j0 -x i0 ||×dv s .w js is the influence weight of the sth neighboring material point in the neighborhood corresponding to the jth neighboring material point on the jth neighboring material point, v sis the volume of the sth neighborhood material point, dv s Indicates v s Integrate, g(j) is the number of neighboring material points corresponding to the jth neighboring material point.

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

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

[0071] In the embodiment of the present invention, an adaptive dynamic relaxation method can be used to convert the near-field dynamics equilibrium equation into an equilibrium equation in the form of an ordinary differential equation by setting virtual damping and virtual mass, so as to solve the displacement of the material point.

[0072] Further, in S400, the displacement x of the material point at the uth time step is u The following conditions must be met:

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

[0074] Among them, x u-1 is the displacement of the u-1th time step, △t is the time step length, that is, the step length between two adjacent time steps, v u+(1 / 2) is the velocity at the u+(1 / 2)th time step, v u+(1 / 2) =((2-dc u ×△t)×v u-(1 / 2) +2×△t×DR -1 ×F u ) / (2+dc u ×△t), dc u is the virtual damping matrix at the uth time step, F u is the internal force on the material point at the uth time step, DR is the virtual diagonal density matrix, v u-(1 / 2) is the velocity at the u-(1 / 2)th time step. In an exemplary embodiment, Δt may be set to 1.

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

[0076] In S500, the stress tensor S at material point i isi The following conditions must be met:

[0077] S i =D:E i ;

[0078] Among them, E i is the strain tensor of 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, represents the tensor product, R is the shape tensor of the material point i after deformation, I is the identity matrix, D is the stiffness matrix, : represents the double dot product, and · represents the dot product.

[0079] In an embodiment of the present invention, the stress tensor may include the principal stresses acting on the x-axis, y-axis and z-axis and the shear stresses acting on the xy, yz and zx planes. In this way, the stress of the material point can be obtained by calculating the stress tensor of the material point, that is, the principal stresses and the shear stresses.

[0080] S600, based on the stress of the material point at the current time step, the Mohr-Coulomb failure criterion is used to judge the fracture condition of the virtual bond, and the local damage value of the material point at the current time step is obtained; if the current time step is equal to the preset time step threshold, exit the current control program, otherwise, execute S300.

[0081] Further, in an embodiment of the present invention, the fracture condition of the virtual bond 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 using the Mohr-Coulomb failure criterion, which can specifically include the following steps:

[0082] S601, obtain the virtual key K corresponding to any material point i j The stress tensor σ at the current time step b ij =(S i +S ij ) / 2,S i is the stress tensor at material point i, S ij is the stress tensor of the jth neighboring material point of material point i.

[0083] S602, obtain K j The normal stress σ Nij =σ b ij n, and obtain K j Shear stress U ij =σ b ij -σ N ij n, 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 K j The corresponding shear stress, τ ij K j Shear strength, τ ij =σ N ij ×tan(φ ij )+c ij ,φ ij K j The mean value of the internal friction angle of the two corresponding material points, c ij K j The mean value of the cohesive force of the two corresponding material points.

[0086] In the embodiment of the present invention, the bond constant is used to indicate the integrity of the bond. If the bond constant is 1, it indicates that the bond is in a broken state, at which point there is no longer an interaction relationship between the two interacting material points, that is, in the embodiment of the present invention, when the shear stress of the bond exceeds the shear strength or the normal stress exceeds 0, the bond breaks, and 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 is the volume v of the jth neighborhood material point j Perform integration.

[0088] S605, if P i = 0, it is determined that the material point i is in an undamaged state, that is, all virtual bonds are not destroyed. If P i =1, it is determined that the material point i is in a completely destroyed state, that is, all virtual bonds are destroyed. If 0<P i <1, it is determined that the material point i is in a state of local destruction, that is, some virtual bonds are destroyed.

[0089] In an embodiment of the present invention, the preset time step threshold may 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 pre-test the influence of different time step thresholds on the stability of the results. For example, the time step thresholds may be set to 50, 200, 500, 1000, and 2000, respectively, and then the maximum displacement obtained by simulation is compared. If the maximum displacement is always changing in the case of 50 to 500, and the displacements obtained by 1000 and 2000 are equal, it is proved that the time step threshold can be set to 1000, which will not consume more model running time and can ensure the stability of the model. In an illustrative embodiment, the time step threshold may be set to 1000.

[0090] In the embodiment of the present invention, the displacement of the material point corresponds to the settlement of the soil, the velocity of the material point corresponds to the settlement velocity of the soil, and the local damage of the material point corresponds to the local damage of the soil, that is, ground fissures. The numerical simulation method of surface deformation evolution based on near-field dynamics provided by the present invention can reveal the physical mechanism of ground settlement disasters turning 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, if Figure 3 As shown in the figure, the length of the geological body model is 2000m, the height is 500m, the thickness is 100m, and the bedrock burial depth is 100m. It is divided into 400, 100, and 20 material points along the length, height, and thickness directions respectively, and the virtual boundary layer is divided into 3 material points with a material point spacing of 5m. The near field range, i.e., the neighborhood, is 3.015 times the material point spacing.

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

[0096] Table 1:

[0097]

[0098] (2) Neighborhood point retrieval:

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

[0100] (3) Boundary conditions:

[0101] The left and right virtual boundary layers of the model are subjected to horizontal fixed constraints, the front and rear virtual boundary layers are subjected to longitudinal fixed constraints, and the lower virtual boundary layer is subjected to horizontal fixed constraints, longitudinal fixed constraints, and vertical fixed constraints. In this embodiment, since the rock is very hard and usually does not cause displacement due to groundwater extraction, the rock area can be treated as a boundary, that is, displacement constraints are directly applied, that is, horizontal fixed constraints, longitudinal fixed constraints, 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 100m. The groundwater level change is used as the power source of the model to solve the force of the material point 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 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 of each virtual bond is calculated. In this embodiment, the time step threshold is set to 1000 steps.

[0104] (5) Damage judgment:

[0105] The fracture of the virtual bond of the material point is judged by the 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, and the local damage value of each material point is obtained by integration.

[0106] (6) Result analysis:

[0107] After the calculation, the evolution process of differential settlement disasters caused by groundwater exploitation under bedrock structural conditions into ground fissures was obtained, such as Figure 4 It can be seen that under the action of groundwater extraction, the soil is compressed, the soil above the bedrock gradually breaks, and finally forms ground fissures.

[0108] In order to verify the effectiveness of the proposed method and device for numerical simulation of surface deformation based on peridynamics, a finite element-interface element model was constructed for comparison. The comparison of the maximum values ​​of horizontal and vertical deformation is shown in Table 2, and the comparison of spatial distribution is shown in Table 2. Figure 5 and Figure 6 shown.

[0109] Table 2:

[0110]

[0111] Embodiment 2:

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

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

[0114] In this embodiment, if Figure 7 As shown in the figure, the length of the geological model is 2000m, the height is 500m, the thickness is 100m, the fault structure is a normal fault structure, and the dip angle is 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 material point spacing is 5m, and the near field range is 3.015 times the material point spacing.

[0115] Initialize each material parameter matrix. The material parameter matrices of this embodiment are the same as those of the first embodiment.

[0116] (2) Neighborhood point retrieval:

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

[0118] (3) Boundary conditions:

[0119] Horizontal fixed constraints are applied to the left and right virtual boundary layers of the model, longitudinal fixed constraints are applied to the front and rear virtual boundary layers, and horizontal fixed constraints, longitudinal fixed constraints and vertical fixed constraints are applied to the lower virtual boundary layer.

[0120] (4) Time domain integration:

[0121] In this implementation case, the groundwater level change of the aquifer is 100m. The groundwater level change is used as the power source of the model to solve the force of the material point 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 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 of each virtual bond is calculated. In this embodiment, the time step threshold is set to 1000 steps.

[0122] (5) Damage judgment:

[0123] The fracture of the virtual bond of the material point is judged by the 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, and the local damage value of each material point is obtained by integration.

[0124] (6) Result analysis:

[0125] After the calculation, the evolution process of differential settlement disasters caused by groundwater exploitation under bedrock structural conditions into ground fissures was obtained, such as 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 close to the fault zone gradually breaks, eventually forming ground fissures.

[0126] In order to verify the effectiveness of the proposed method and device for numerical simulation of surface deformation based on peridynamics, a finite element-interface element model was constructed for comparison. The comparison of the maximum values ​​of horizontal and vertical deformation is shown in Table 3, and the spatial distribution is shown in Table 3. Fig. 9 and Fig.10 shown.

[0127] Table 3:

[0128]

[0129]

[0130] Embodiment three:

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

[0132] (1) Model discretization:

[0133] In this embodiment, if Fig.11 As shown in the figure, the length of the geological model is 2000m, the height is 500m, the thickness is 100m, and the outcrop dip angles are 75° and 105° respectively. It is divided into 400, 100, and 20 material points along the length, height, and thickness directions respectively, and 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 each material parameter matrix. The material parameter matrices of this embodiment are the same as those of the first embodiment.

[0135] (2) Neighborhood point retrieval:

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

[0137] (3) Boundary conditions:

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

[0139] (4) Time domain integration:

[0140] In this implementation case, the groundwater level change of the aquifer is 100m. The groundwater level change is used as the power source of the model to solve the force of the material point 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 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 of each virtual bond is calculated. In this embodiment, the time step threshold is set to 1000 steps.

[0141] (5) Damage judgment:

[0142] The fracture of the virtual bond of the material point is judged by the 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, and the local damage value of each material point is obtained by integration.

[0143] (6) Result analysis:

[0144] After the calculation, the evolution process of differential settlement disasters caused by groundwater exploitation under bedrock structural conditions into ground fissures was obtained, such as Fig.12 As shown in Figure 2, it can be seen that under the action of groundwater extraction, the soil is compressed, the soil close to the bedrock gradually breaks, and finally forms ground fissures.

[0145] In order to verify the effectiveness of the proposed surface deformation numerical simulation method and system based on peridynamics, a finite element-interface element model was constructed for comparison. The comparison of the maximum values ​​of horizontal and vertical deformation is shown in Table 4, and the comparison of spatial distribution is shown in Table 4. Fig.13 and Fig.14 shown.

[0146] Table 4:

[0147]

[0148] The above three embodiments have achieved simulation effects close to those of the traditional finite element model. It can be seen that the numerical simulation method of surface deformation evolution based on near-field dynamics provided by the present invention can effectively simulate the process of differential ground subsidence disasters caused by groundwater exploitation under various geological structural conditions such as bedrock, fractures and outcrops turning into ground fissures, while avoiding the defect of the finite element-interface element model that it is necessary to know in advance whether cracks exist and their positions and sizes.

[0149] Based on the same inventive concept, an embodiment of the present invention provides a numerical simulation device for surface deformation evolution based on peridynamics, the device comprising:

[0150] The initial geological body model acquisition module is used to discretize the 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 on the outside of 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 includes the number, coordinates and material parameters of the material point.

[0151] The material point neighborhood acquisition module is used to apply corresponding boundary conditions on the virtual boundary layer and acquire the neighborhood of each material point in the initial geological body model, wherein the distance between any material point and any neighboring material point in its corresponding neighborhood is less than a set distance d0.

[0152] The material point force acquisition module is used to obtain the internal force received 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 received by the material point at the current time step.

[0153] The material point displacement acquisition module is used to obtain the displacement of the material point in the current time step based on the force on the material point in the current time step using the adaptive dynamic relaxation method.

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

[0155] The material point damage judgment module is used to judge the fracture of the virtual bond based on the stress of the material point at the current time step and the Mohr-Coulomb failure criterion, and obtain the local damage value of the material point at the current time step. Figure 1 Therefore, for the functions that can be realized by each functional module of the device, please refer to Figure 1 The description of the illustrated embodiment is not repeated in detail.

[0156] An embodiment of the present invention also provides an electronic device, comprising: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, and the instructions are configured to execute the method described in the embodiment of the present invention.

[0157] The embodiment of the present invention further provides a non-transitory computer-readable storage medium storing computer-executable instructions, wherein the computer instructions are used to execute the method described in the embodiment of the present invention.

[0158] It should be understood that the various forms of processes shown above can be used to reorder, add or delete steps. For example, the steps described in the present invention can be executed in parallel, sequentially or in different orders, as long as the desired results of the technical solution disclosed in the present invention can be achieved, and this document does not limit this.

[0159] The above specific implementations do not constitute a limitation on the protection scope of the present invention. It should be understood by those skilled in the art that various modifications, combinations, sub-combinations and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for numerical simulation of surface deformation evolution based on peridynamics, characterized in that: The method comprises the following steps: S100, discretizing the geological body model to be calculated into n1 material points, obtaining corresponding material point information, 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 includes the number, coordinates and material parameters of the material point; S200, applying corresponding boundary conditions on the virtual boundary layer, and obtaining the neighborhood of each material point in the initial geological body model, wherein the distance between any material point and any neighboring material point in its corresponding neighborhood is less than a set distance d0; S300, obtaining the internal force received 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 received by the material point at the current time step; S400, based on the force of the material point at the current time step, using an adaptive dynamic relaxation method to obtain the displacement of the material point at the current time step; S500, based on the displacement of the material point at the current time step, obtaining the stress of the material point at the current time step; S600, based on the stress of the material point at the current time step, use the Mohr-Coulomb failure criterion to judge the fracture condition of the virtual bond, and 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, exit the current control program, otherwise, execute S300.

2. The method according to claim 1, characterized in that: d0=f×Δd, f is a preset value greater than 1, and Δd is the distance between two adjacent material points in the initial geological body model.

3. The method according to claim 1, characterized in that The interaction relationship between any material point and any neighboring material point in the neighborhood corresponding to the material point is represented by a virtual bond; The internal force F acting on any material point i i The following conditions are met: F i =∫ h(i) j＝1 (F ij -F ji )×dv j ; Among them, F ij is the force exerted by the jth neighboring material point in the neighborhood corresponding to material point i on material point i, F ji is the force exerted by material point i on the jth neighboring material point in its neighborhood, i ranges from 1 to (n1+n2), j ranges from 1 to h(i), h(i) is the number of neighboring material points in the neighborhood corresponding to material point i, v j is the volume of the jth neighborhood material point, dv j Indicates v j Perform integration; Among them, F ji The following conditions must be met: F ji =((-3)×p i / wv i )×w ji ×||x j0 -x i0 ||+(15×SM i ) / wv i ×w ji ×e ij d ; w ji is the influence weight of material point i on the jth neighboring material point, |||| represents the norm, x i0 is the initial coordinate of material point i, i.e., the coordinate in the initial geological model, x j0 is the initial coordinate of the jth neighborhood material point, SM i is the shear modulus of material point i, e ij d is the elongation state e of the virtual bond between material point i and the jth neighboring material point ij The deviation component, wv i is the weighted volume of material point i, p i is the peridynamic pressure at material point i, p i =-k i ×θ i +γ i ×pf i , k i is the bulk modulus of material point i, γ i is the fluid pressure coefficient at material point i, pf i is the change in pore water pressure at material point i, θ i is the volume expansion rate of material point i; F ij The following conditions must be met: F ij =((-3)×p j / wv j )×w ij ×||x j0 -x i0 ||+(15×SM j ) / wv j ×w ij ×e ij d ; Among them, w ij is the influence weight of the jth neighboring material point on material point i, SM j is the shear modulus of the jth neighboring material point, wv j is the weighted volume of the jth neighborhood material point, p j is the peridynamic pressure of the jth neighboring material point, p j =-k j ×θ j +γ j ×pf j , k j is the bulk modulus of the jth neighborhood material point, γ j is the fluid pressure coefficient of the jth neighboring material point, pf j is the change in pore water pressure of the jth neighboring material point, θ j is the volume expansion rate of the jth neighborhood material point.

4. The method according to claim 3, characterized in that 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 is the displacement of material point i, x jc is the displacement of the jth neighborhood material point; wv i =∫1 h(i) w ij ×||x j0 -x i0 ||×||x j0 -x i0 ||×dv j ; θ i =(3 / wv i )×∫1 h(i) w ij ×||x j0 -x i0 ||·e ij ×dv j 。 5. The method according to claim 1, characterized in that In S400, the displacement x of the material point at the uth time step is u The following conditions must be met: x u =x u-1 +v u+(1 / 2) ×△t; Among them, x u-1 is the displacement of the u-1th time step, △t is the time step, v u+(1 / 2) is the velocity at the u+(1 / 2)th time step, v u +(1 / 2) =((2-dc u ×△t)×v u-(1 / 2) +2×△t×DR -1 ×F u ) / (2+dc u ×△t), dc u is the virtual damping matrix at the uth time step, F u is the internal force on the material point at the uth time step, DR is the virtual diagonal density matrix, v u-(1 / 2) is the velocity at the u-(1 / 2)th time step.

6. The method according to claim 3, characterized in that In S500, the stress tensor S at any material point i is i The following conditions must be met: S i =D:E i ; E i is the strain tensor of 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, represents the tensor product, R is the shape tensor of the material point i after deformation, I is the identity matrix, D is the stiffness matrix, : represents the double dot product, and · represents the dot product.

7. The method according to claim 5, characterized in that The method of determining the fracture of the virtual bond based on the stress of the obtained material point and using the Mohr-Coulomb failure criterion to obtain the local damage value of the material point includes the following steps: S601, obtain the virtual key K corresponding to the material point i j The stress tensor σ b ij =(S i +S ij ) / 2,S i is the stress tensor at material point i, S ij is the stress tensor of the jth neighboring material point of material point i, and the value of j ranges from 1 to h(i); S602, obtain K j The normal stress σ N ij =σ b ij n, and obtain K j Shear stress U ij =σ b ij -σ N ij n, n is the direction vector; S603, if ||U ij ||≥τ ij Or N ij ≥0, set K j The corresponding bond constant g ij =1, if ||U ij ||<τ ij And N ij <0, set g ij =0,U ij is the shear stress corresponding to virtual bond j, τ ij K j Shear strength, τ ij =σ N ij ×tan(φ ij )+c ij ,φ ij K j The mean value of the internal friction angle of the two corresponding material points, c ij K j The mean value of the cohesive force of the two corresponding material points; S604, obtaining the local damage value P of material point i i =1-(∫1 h(i) g ij ×dv j / ∫1 h(i) dv j ); S605, if P i = 0, it is determined that the material point i is in an undamaged state, if P i =1, it is determined that the material point i is in a completely destroyed state, if 0<P i <1, it is determined that the material point i is in a state of local destruction.

8. The method according to claim 1, characterized in that The virtual boundary layers 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; Among them, applying corresponding boundary conditions on the virtual boundary layer includes: applying horizontal fixed constraints to the left virtual boundary layer and the right virtual boundary layer, applying longitudinal fixed constraints to the front virtual boundary layer and the rear virtual boundary layer, and applying horizontal fixed constraints, longitudinal fixed constraints and vertical fixed constraints to the lower virtual boundary layer.

9. A numerical simulation device for surface deformation evolution based on peridynamics, characterized in that: The device comprises: The initial geological body model acquisition module is used to discretize the geological body model to be calculated into n1 material points, obtain the corresponding material point information, 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 includes the number, coordinates and material parameters of the material point; A material point neighborhood acquisition module, used to apply corresponding boundary conditions on the virtual boundary layer and acquire the neighborhood of each material point in the initial geological body model, wherein the distance between any material point and any neighboring material point in its corresponding neighborhood is less than a set distance d0; A material point force acquisition module is used to obtain the internal force of the material point in the current time step based on the change of the pore water pressure of the material point in the current time step, as the force of the material point in the current time step; The material point displacement acquisition module is used to obtain the displacement of the material point at the current time step by using the adaptive dynamic relaxation method based on the force of the material point at the current time step; Material point stress acquisition module, used to obtain the stress of the material point at the current time step based on the displacement of the material point at the current time step; The material point damage judgment module is used to judge the fracture of the virtual bond based on the stress of the material point at the current time step and the Mohr-Coulomb failure criterion to obtain the local damage value of the material point at the current time step.

Images