Approximate simulation method of small deformation mode of FLAC3D based on displacement and stress correction

Abstract

The application discloses a kind of FLAC3D small deformation mode approximate simulation method based on displacement and stress correction, numerical simulation field.Based on the small deformation mode of FLAC3D software, the target working condition is simulated initially;Extract displacement nephogram, judge grid and node embedding pair, determine displacement limited area according to the displacement that does not comply with actual area;Extract the pointer of the node and the grid in the area, write in list, process the grid and node pointer that mutually embed at boundary, form grid and node embedding pair list;Combined with the displacement interval of actual object, write displacement and stress correction function;The target working condition is simulated twice and the displacement and stress in the small deformation mode simulation process are corrected;Repeat the above steps until there is no mutual embedding after grid point coordinates are superposed with displacement parameters.The method can reduce modeling requirements, adapt to more complex numerical models, and better reflect the stress, displacement, plastic zone and other distribution of actual engineering.

Description

Technical Field

[0001] This invention belongs to the field of numerical simulation, specifically relating to an approximate simulation method for small deformation modes in FLAC3D based on displacement and stress correction. Background Technology

[0002] Current scientific research often employs methods such as theoretical analysis, laboratory experiments, numerical simulation, and field practice. Among these, laboratory experiments require specialized equipment and incur certain experimental costs, while field practice is large-scale and requires coordination and collaboration from various parties. Numerical simulation, on the other hand, has advantages such as lower research costs, easier result acquisition, the ability to perform multi-factor and multi-level analysis, and the ability to perform multi-condition analysis at the engineering scale, making it a favorite among many scholars.

[0003] FLAC3D, a commonly used numerical simulation software in geotechnical engineering, is based on the finite difference method and boasts advantages such as accurate and reasonable calculation results, applicability to unstable processes, small memory footprint, and fast computation speed. The FLAC3D calculation process can be divided into two modes: large deformation and small deformation. In the large deformation mode, the mesh node coordinates change in real time to reach an equilibrium state, resulting in more accurate calculation results and wider applicability, but the computational cost is higher, and mesh deformation may cause program errors. In the small deformation mode, the node coordinates remain unchanged, resulting in fast computation speed and fewer program errors, but it is not suitable for simulating large deformations. Summary of the Invention

[0004] To address the shortcomings of existing technologies, a method for approximate simulation of small deformation modes in FLAC3D based on displacement and stress correction is provided. This method is simple to implement and easy to use. By simulating the working conditions and correcting the displacement and stress during the small deformation mode simulation process, the grid point coordinates are superimposed with displacement parameters without interlocking regions, effectively avoiding errors caused by grid deformation. The simulation results for working conditions with large deformations are more realistic than those for general small deformation calculations.

[0005] To achieve the above technical objectives, this invention discloses an approximate simulation method for small deformation modes in FLAC3D based on displacement and stress correction, specifically including the following steps:

[0006] s1. Establish an initial numerical model in which the grid points on the surface of the mining face and the roof and floor of the tunnel overlap on the horizontal plane. The numerical model consists of grids and grid points. All grids and grid points have their own ID numbers. Based on the small deformation mode of FLAC3D software, perform a preliminary simulation of the displacement changes of the rock strata and surrounding rock during mining and tunneling. Record the total number of simulation steps as S.

[0007] s2. Extract the coordinates of grid points before and after displacement from the preliminary simulation results. The coordinates before displacement are the initial coordinates, and the coordinates after small deformation displacement are the final coordinates. The initial coordinates of the grid points are affected by the displacement parameters generated during mining. The final coordinates of all grid points on the lower surface of the roof and the bottom rock strata in the numerical model are calculated using the formula A'=A+D, where A' represents the final coordinates of the grid point, A represents the initial coordinates of the same grid point, and D represents the displacement of the same grid point. The final coordinates of all grid points on the lower surface of the roof and the bottom rock strata are plotted using the above formula. The final coordinates of the grid points on the bottom rock strata are located above the lower surface of the roof and are denoted as the bottom embedded region. Similarly, the final coordinates of the grid points on the upper surface of the bottom and all grid points on the roof rock strata are calculated using the formula A'=A+D. All grid points on the lower surface of the roof and the bottom rock strata are plotted based on the final coordinates.

[0008] If the lower surface of the top plate and the upper surface of the bottom plate meet during a small deformation simulation and are not in equilibrium, they will continue to move along their respective displacement directions to form an embedding. Therefore, the area where the final coordinates of all grid points of the top plate rock layer are below the upper surface of the bottom plate is recorded as the top plate embedding area. Using the final coordinates of all grid points of the top and bottom plate rock layers, the lower surface of the top plate and the upper surface of the bottom plate are drawn. The area where the lower surface of the top plate and the upper surface of the bottom plate intersect is the area where the top and bottom plates are mutually embedded. The coordinates of each grid point in the area where the top and bottom plates are mutually embedded are obtained. The three-axis data of the coordinates of each grid point in the area where the top and bottom plates are mutually embedded are sorted by size to obtain the maximum and minimum values ​​of the three-axis coordinates. The area between the maximum and minimum values ​​of the three-axis coordinates is determined as the displacement-limited area.

[0009] s3. Based on the displacement-limited area obtained in step s2, extract the initial coordinates of the surface grid and grid points of the mining face and the top and bottom plates of the tunnel within the displacement-limited area, including the initial coordinates of the grid on the lower surface of the top plate, the initial coordinates of the grid on the upper surface of the bottom plate, the initial coordinates of the grid points on the lower surface of the top plate, and the initial coordinates of the grid points on the upper surface of the bottom plate.

[0010] Within the displacement-limited area, the coordinates of the grid on the lower surface of the top plate are used as elements to generate a top plate grid list in ascending order of the corresponding grid ID in the numerical model. The elements in the top plate grid list are sorted, and the coordinates of the bottom plate grids that coincide with the horizontal plane of each top plate grid are sorted and used as elements to generate a bottom plate grid list. The coordinates of all grid points on the lower surface of the top plate within the displacement-limited area are sorted and used as elements to generate a top plate grid point list. The elements in the top plate grid point list are sorted, and the coordinates of the bottom plate grid points that coincide with the horizontal plane of each top plate grid point are sorted and used as elements to generate a bottom plate grid point list.

[0011] Elements with the same ID in the top and bottom grid lists are treated as grid embedding pairs, and elements with the same ID in the top and bottom grid point lists are treated as grid point embedding pairs. The value of the ID is the number of list elements after the data is added.

[0012] s4. Based on mesh embedding pairs and mesh point embedding pairs, a small deformation secondary simulation is performed on the initial numerical model. During the simulation, the numerical model information of the small deformation secondary simulation is extracted every 1-1 / 2 s steps, including the displacement vectors of the mesh points represented by the elements in the top and bottom plate mesh point lists. Each displacement vector is added as an element to the real-time top plate displacement list. According to the recording order of the top plate mesh points corresponding to each element in the real-time top plate displacement list, the order of the corresponding bottom plate mesh points is obtained according to the mesh point embedding pair relationship. The displacement parameters of each bottom plate mesh point are extracted according to the bottom plate mesh point order and added to the real-time bottom plate displacement list. The real-time surrounding rock deformation list is obtained by subtracting the corresponding element in the real-time bottom plate displacement list from each element in the real-time top plate displacement list. The positions of elements whose deformation components in the surrounding rock deformation list are greater than the dimension of this component in the excavated area are mapped to the real-time top and bottom plate displacement list. This yields the real-time displacement limit area in the process of calculating the unbalanced force to the target value using the numerical model of small deformation secondary simulation. Displacement correction is achieved by fixing the velocity of the large displacement side grid points in the top and bottom plates within the real-time displacement limit area to 0. The displacement correction in the real-time displacement limit area is removed using FLAC3D software commands, and an equivalent reaction force is applied to the grid points in the area where the displacement correction is removed. The reaction force of the grid points in the original real-time displacement limit area is extracted. Combining the embedding relationship between the top plate grid points and the bottom plate grid points, pressure is applied to the corresponding grid points in the top plate when extracting the bottom plate grid points, and pressure is applied to the corresponding grid points in the bottom plate when extracting the corresponding grid points in the top plate.

[0013] After the pressure is applied, the displacement on the stress extraction side is fixed to 0. The mesh stress of the shallow layer of the self-excavated area on the displacement-limited side is extracted using FLAC3D software. The shallow layer is the first stratum outward from the top and bottom plates of the excavated area, and the thickness is set according to actual needs. An approximate calculation formula is used. Transformation, combined with mesh embedding, corrects the mesh stress of the first layer of the bottom plate corresponding to the top plate or the bottom plate corresponding to the top plate from the excavated area outward, and generates displacement and stress correction functions between rock layers or surrounding rocks.

[0014] Where, σ i ’ The stress in the shallow mesh on the side to be corrected is given by: σ0 is the stress in the shallow mesh on the side to be corrected as calculated in the simulation before correction; E2 is the Young's modulus of the side to be corrected; E1 is the Young's modulus of the shallow mesh on the displacement-constrained side; μ1 is the Poisson's ratio of the shallow mesh on the displacement-constrained side; and μ2 is the Poisson's ratio of the shallow mesh on the side to be corrected. j σ j To simulate the stress in the shallow mesh under displacement constraints after a certain number of steps, σ j0 To simulate the stress in the shallow layer of the side-shallow mesh within a certain number of steps of pre-displacement, ε i The strain of the shallow mesh on the side to be corrected;

[0015] s5. Repeat steps s2-s4 above until the grid point coordinates of the rock strata and surrounding rock are superimposed with displacement parameters and there is no mutual embedding.

[0016] Furthermore, the displacement-limited region is larger than the region where displacement parameters are inter-embedded after the coordinates of the small deformation simulation mesh are superimposed.

[0017] Furthermore, during the implementation of the method, displacement correction includes three processes: traversal, judgment, and restriction, while stress correction includes three processes: reaction force extraction, pressure application, and stress correction.

[0018] Furthermore, the displacement correction process traverses the entire list of grid points and all elements of the grid list in the numerical model; the judgment process determines the range by comparing the sum of the original coordinates and real-time displacement of the grid point with the actual position limit; the constraint process is achieved by fixing the node velocity in one or more directions to 0, and deleting grid points in the list whose displacement has been constrained.

[0019] Furthermore, the method for obtaining the stress correction function is as follows: Reaction extraction involves using the "zone gridpoint force-reaction" command to remove displacement constraints and apply equivalent reaction forces, then using a loop to iterate through the grid points within the displacement-limited region to extract the reaction forces at the nodes in the original restricted region. Pressure application involves calculating the approximate stress in the vicinity of each grid point using the reaction forces at the grid points on the displacement-fixed side and the area of ​​the adjacent region. Using the "zone face apply stress" command, the grid points on the non-displacement-fixed side are obtained through grid point embedding relationships. Stress is applied to the grid points in contact with the non-displacement-fixed side. After pressure application, the displacement on the stress extraction side is reset to 0. Stress correction involves extracting the stress from the shallow grid of the top plate and applying an approximate calculation formula. The transformation, combined with mesh embedding, is a correction function for the mesh stress of the first layer of the bottom plate from the excavated region outwards, where σ i ’ Shallow stress in the base plate, σ0 is the uncorrected shallow stress in the base plate, E2 is the Young's modulus of the base plate, E1 is the Young's modulus of the top plate, μ1 is the Poisson's ratio of the top plate, μ2 is the Poisson's ratio of the base plate, σ j σ j To simulate the stress in the shallow top plate after a certain number of steps, σ j0 To simulate the stress in the shallow layer of the top plate before a certain number of steps, ε i This represents the shallow strain of the base plate.

[0020] Furthermore, after the numerical model of the small deformation secondary simulation is self-equilibrium, the "model step" command is used for solving. According to the parameters corresponding to the preset number of numerical simulation steps, the displacement correction function and stress correction function must be called after each line of solution command. When the numerical model is not balanced, the displacement correction function and stress correction function are called once after each preset running time.

[0021] A computer device includes a processor and a memory, the processor being electrically connected to the memory, the memory being used to store instructions and data, and the processor being used to execute an approximate simulation method for small deformation modes in FLAC3D based on displacement and stress correction.

[0022] Beneficial Effects: This invention provides an approximate simulation method for small deformation modes in FLAC3D based on displacement and stress correction. First, based on the FLAC3D small deformation mode, a preliminary simulation of the target working condition is performed, displacement contour maps are extracted, mesh and node embedding pairs are determined, and displacement-limited regions are identified based on areas where displacement does not conform to reality. Second, pointers to nodes and meshes within these regions are extracted and written into a list. Pointers to mutually embedded meshes and nodes at the boundaries are processed to form a list of mesh and node embedding pairs. Then, combining with real-world conditions, displacement and stress correction functions are written to perform a secondary simulation of the target working condition and correct the displacement and stress during the small deformation mode simulation process. Finally, the above steps are repeated until there are essentially no areas where displacement does not conform to reality. This method can avoid program errors caused by mesh deformation, and the simulation results for working conditions with large deformations are more realistic than general small deformation calculations. Compared with existing technologies, it has the following advantages:

[0023] 1) This invention expands the applicability of small deformation modes and provides more accurate simulation results.

[0024] 2) This invention reduces the modeling requirements for simulating large deformation conditions.

[0025] 3) This invention avoids program errors caused by mesh deformation during simulation. Attached Figure Description

[0026] Figure 1 This is a flowchart of the FLAC3D small deformation mode approximation simulation method based on displacement and stress correction according to the present invention.

[0027] Figure 2 This is a displacement contour map of the z-direction in the small deformation mode in an embodiment of the present invention.

[0028] Figure 3 This is a stress cloud diagram in the z-direction after displacement and stress correction in an embodiment of the present invention.

[0029] Figure 4 This is a displacement cloud diagram in the z-direction after displacement and stress correction in an embodiment of the present invention.

[0030] Figure 5 This is a diagram showing the distribution of the plastic zone after displacement and stress correction in an embodiment of the present invention. Detailed Implementation

[0031] The embodiments of the present invention will be further described below with reference to the accompanying drawings:

[0032] This invention discloses an approximate simulation method for small deformation modes in FLAC3D based on displacement and stress correction:

[0033] s1. Based on the small deformation mode of FLAC3D software, a preliminary simulation of the target working condition is performed;

[0034] s2. Based on the preliminary simulation results of the target working condition, construct the relationship between the top grid and grid point pairs: extract the displacement cloud map, determine the grid and node embedding pairs, and determine the displacement-limited area based on the mutual embedding area after superimposing displacement parameters on the grid point coordinates. The displacement-limited area is larger than the mutual embedding area after superimposing displacement parameters on the grid point coordinates.

[0035] s3. Based on the displacement-limited region and mesh / node embedding pairs from step s2, extract the pointers of nodes and meshes within the region, write them into a list, process the mesh and node pointers that are mutually embedded at the boundary, and form a list of mesh / node embedding pairs.

[0036] s4. Based on the pointer list in step s3, and combined with real-world scenarios, write a displacement correction function that includes traversal, judgment, and restriction, and a stress correction function that includes reaction force extraction, pressure application, and stress correction.

[0037] s5. Based on the displacement and stress correction function in step s4, perform a second simulation of the target working condition and after the model self-equilibrates during the simulation of the small deformation mode, call the displacement and stress correction function after each line of solution command. When the model is not balanced, call the displacement and stress correction function once after calculating 1 to 1 / 2S (S is the total number of steps in the initial simulation) steps to make corrections.

[0038] s6. Repeat steps s2-s6 until there are no inter-intercalated areas after the grid point coordinates are superimposed with the displacement parameters.

[0039] like Figure 1 As shown, a method for approximate simulation of small deformation modes in FLAC3D based on displacement and stress correction includes the following steps:

[0040] s1. Establish a numerical model where the grid points on the roof and floor surfaces of the longwall face and tunneling roadway coincide on the horizontal plane. Based on the FLAC3D small deformation mode, perform a preliminary simulation of the displacement changes of the rock strata and surrounding rock during longwall face mining and tunneling. Record the total number of simulation steps as S; the z-direction displacement cloud map of the small deformation mode is shown below. Figure 2 As shown;

[0041] s2. Extract the initial coordinates of the grid points and the displacement parameters of the grid points caused by mining from the preliminary simulation results. Calculate the final coordinates of all grid points on the lower surface of the roof and the bottom strata using the formula A'=A+D, where A' is the final coordinate of the grid point, A is the initial coordinate of the same grid point, and D is the displacement of the same grid point. Draw all grid points on the lower surface of the roof and the bottom strata based on the final coordinates. The final coordinates of the bottom strata grid points are above the lower surface of the roof, denoted as the bottom embedded region. Calculate the final coordinates of the upper surface of the bottom strata and all grid points on the roof strata using the formula A'=A+D. Draw all grid points on the lower surface of the roof and the bottom strata based on the final coordinates. Due to the small deformation simulation of the roof and bottom... When the surfaces of the plates meet, if they are not in equilibrium, they will continue to move along their respective displacement directions to form an embedding. Therefore, the area where the final coordinates of the top plate strata are below the top surface of the bottom plate is recorded as the embedding region of the top plate. The lower surface of the top plate and the upper surface of the bottom plate are drawn using coordinates. The area where the lower surface of the top plate and the upper surface of the bottom plate intersect is the embedding region between the top and bottom plates. The coordinates of each point in the embedding region between the top and bottom plates are obtained. The three-axis data of the coordinates of each point in the embedding region between the top and bottom plates are sorted by size to obtain the maximum and minimum values ​​of the three-axis coordinates. The area between the maximum and minimum values ​​of the three-axis coordinates is determined as the displacement-limited region. The displacement-limited region is larger than the embedding region after the displacement parameters are superimposed on the coordinates of the small deformation simulation mesh.

[0042] s3. Based on the displacement-limited area obtained in step s2, extract the initial coordinates of the surface grids and grid points of the mining face and the top and bottom plates of the tunnel within the displacement-limited area. Add all grid coordinates of the lower surface of the top plate within the displacement-limited area to the top plate grid list. Set the order in which the grid coordinates of the lower surface of the top plate are added to the top plate grid list to the element order of the top plate grid point list. Add the bottom plate grid coordinates that coincide with the horizontal plane of each top plate grid to the bottom plate grid list according to the element order of the top plate grid list. Add all grid point coordinates of the lower surface of the top plate within the area to the top plate grid point list. Add the bottom plate grid point coordinates that coincide with the horizontal plane of each top plate grid point to the bottom plate grid point list according to the element order of the top plate grid point list. Record the grid coordinate elements with the same position and number in the top plate grid list and the bottom plate grid list as grid embedding pairs. Record the grid points with the same ID number in the top plate grid point list and the bottom plate grid point list as grid point embedding pairs.

[0043] s4. Based on the mesh and mesh point embedding pairs from step s3, perform a secondary step-by-step simulation with small deformations on the original model from step s1. During the secondary step-by-step simulation, extract the displacement parameters of the mesh points represented by the elements in the top plate mesh point list of the model in the secondary simulation of small deformations every 1-1 / 2 seconds. Construct a real-time top plate displacement list for each displacement parameter. The displacement parameters in the real-time top plate displacement list are recorded in the following order: the displacement parameters in the real-time top plate displacement list exist in order 'a' according to their id number (addition time), and a set of displacement parameters (x, y, z) originates from a top plate mesh point 'b', so the top plate mesh points also have an order 'a1'. One top plate mesh point corresponds to one bottom plate mesh point 'b'. When processing the bottom plate mesh points, the order 'a1' of 'b' is used. The algorithm is as follows: Based on the mesh point embedding relationship, the order of the corresponding bottom plate mesh points is obtained. The displacement parameters of each bottom plate mesh point are extracted according to the bottom plate mesh point order. The displacement parameters are added to the real-time bottom plate displacement list. The elements of the real-time top plate displacement list are subtracted from the corresponding elements in the real-time bottom plate displacement list to obtain the real-time surrounding rock deformation list. The position of the element whose deformation component of each element in the real-time surrounding rock deformation list is greater than the dimension of the component in the excavated area is recorded and mapped to the real-time top and bottom plate displacement list. The real-time displacement limit area in the process of calculating the unbalanced force to the target value in the model of small deformation secondary simulation is obtained. The displacement correction is achieved by fixing the velocity of the large displacement side in the top and bottom plates within the real-time displacement limit area to 0. After the numerical model of the secondary simulation is self-balanced, the "model step" command is used for solving. According to the parameters corresponding to the preset number of simulation steps, the displacement and stress correction functions must be called after each line of solution command. When the model is not balanced, the displacement and stress correction functions are called once after each preset running time.

[0044] The FLAC3D software commands are used to remove the displacement restriction of the real-time displacement-limited area and apply equivalent reaction forces. The reaction forces of the nodes in the original real-time displacement-limited area are extracted. Based on the embedding relationship between the top and bottom plate mesh points, pressure is applied to the corresponding top plate mesh points when a bottom plate mesh point is extracted, and vice versa. After applying pressure, the displacement on the stress extraction side is reset to 0. The shallow mesh stress on the displacement-limited side is then extracted. The shallow layer is the first stratum outward from the top and bottom plates of the excavated area, and its thickness is set according to actual needs. Approximate calculation formulas are used. Transformation, combined with mesh embedding, corrects the stress of the shallow mesh of the top plate corresponding to the bottom plate or the bottom plate corresponding to the top plate in the list, and generates displacement and stress correction functions between rock layers or surrounding rocks.

[0045] In FLAC3D software, the displacement correction function consists of three parts: traversal, judgment, and restriction. The stress correction function consists of three parts: reaction force extraction, pressure application, and stress correction. The displacement correction function traverses all elements in the node and mesh pointer list. The judgment part determines the position by comparing the sum of the original coordinates and real-time displacement of the node with the actual position limit. The restriction part is achieved by fixing the node velocity in one or more directions to 0 and deleting node pointers in the list whose displacement has been restricted.

[0046] The method for obtaining the stress correction function is as follows: Reaction force extraction involves removing displacement constraints and applying equivalent reaction forces using the "zone gridpoint force-reaction" command, then iterating through the node pair list using a loop to extract the reaction forces of nodes in the original constrained region. Pressure application involves calculating the approximate stress in the vicinity of each node using the node reaction forces and the area of ​​the adjacent region. Using the "zone face apply stress" command, combined with the embedded pair list, stress is applied to the vicinity of the side in contact with the original node / mesh. After pressure application, the displacement on the stress extraction side is reset to 0. Stress correction involves extracting the stress from the shallow mesh of the top plate and applying an approximate calculation formula. The transformation, combined with the mesh embedding, corrects the stress of the shallow mesh of the bottom plate by a function, where σ i ’ Shallow stress in the base plate, σ0 is the uncorrected shallow stress in the base plate, E2 is the Young's modulus of the base plate, E1 is the Young's modulus of the top plate, μ1 is the Poisson's ratio of the top plate, μ2 is the Poisson's ratio of the base plate, σ j σ j To simulate the stress in the shallow top plate after a certain number of steps, σ j0 To simulate the stress in the shallow layer of the top plate before a certain number of steps, ε i This represents the shallow strain of the base plate.

[0047] Where, σ i ’ The stress in the shallow mesh on the side to be corrected is given by: σ0 is the stress in the shallow mesh on the side to be corrected as calculated in the simulation before correction; E2 is the Young's modulus of the side to be corrected; E1 is the Young's modulus of the shallow mesh on the displacement-constrained side; μ1 is the Poisson's ratio of the shallow mesh on the displacement-constrained side; and μ2 is the Poisson's ratio of the shallow mesh on the side to be corrected. j σ j To simulate the stress in the shallow mesh under displacement constraints after a certain number of steps, σ j0 To simulate the stress in the shallow layer of the side-shallow mesh within a certain number of steps of pre-displacement, ε i The strain of the shallow mesh on the side to be corrected;

[0048] s5. Repeat steps s2-s4 above until the grid point coordinates of the rock strata and surrounding rock are superimposed with displacement parameters and there is no mutual embedding; effectively reflecting the distribution of stress, displacement, plastic zone, etc. in actual engineering, as shown in step 3. Figure 4 and Figure 5 As shown.

Claims

1. A method for approximate simulation of small deformation modes in FLAC3D based on displacement and stress correction, characterized in that, Specifically, the following steps are included: s1. Establish an initial numerical model in which the grid points on the surface of the mining face and the tunnel roof and floor overlap on the horizontal plane. The numerical model consists of grids and grid points. All grids and grid points have their own ID numbers. Based on the small deformation mode of FLAC3D software, perform a preliminary simulation of the displacement changes of the rock strata and surrounding rock during mining and tunnel excavation. Record the total number of simulation steps as S. s2. Extract the coordinates of grid points before and after displacement from the preliminary simulation results. The coordinates before displacement are the initial coordinates, and the coordinates after small deformation displacement are the final coordinates. The initial coordinates of the grid points are affected by the displacement parameters generated by mining. The final coordinates of all grid points on the lower surface of the roof and the bottom rock strata in the numerical model are calculated using the formula A'=A+D. In the formula, A' represents the final coordinates of the grid point, A represents the initial coordinates of the same grid point, and D represents the displacement of the same grid point. The final coordinates of all grid points on the lower surface of the top plate and the bottom plate rock layer are drawn using the above formula. The final coordinates of the grid points of the bottom plate rock layer are above the lower surface of the top plate and are denoted as the bottom plate embedded area. Similarly, the final coordinates of the grid points on the upper surface of the bottom plate and all grid points of the top plate rock strata are calculated using the formula A'=A+D, and all grid points of the lower surface of the top plate and the bottom plate rock strata are drawn based on the final coordinates. If the lower surface of the top plate and the upper surface of the bottom plate meet during a small deformation simulation and are not in equilibrium, they will continue to move along their respective displacement directions to form an embedding. Therefore, the area where the final coordinates of all grid points of the top plate rock layer are below the upper surface of the bottom plate is recorded as the top plate embedding area. Using the final coordinates of all grid points of the top and bottom plate rock layers, the lower surface of the top plate and the upper surface of the bottom plate are drawn. The area where the lower surface of the top plate and the upper surface of the bottom plate intersect is the area where the top and bottom plates are mutually embedded. The coordinates of each grid point in the area where the top and bottom plates are mutually embedded are obtained. The three-axis data of the coordinates of each grid point in the area where the top and bottom plates are mutually embedded are sorted by size to obtain the maximum and minimum values ​​of the three-axis coordinates. The area between the maximum and minimum values ​​of the three-axis coordinates is determined as the displacement-limited area. s3. Based on the displacement-limited area obtained in step s2, extract the initial coordinates of the surface grid and grid points of the mining face and the top and bottom plates of the tunnel within the displacement-limited area, including the initial coordinates of the grid on the lower surface of the top plate, the initial coordinates of the grid on the upper surface of the bottom plate, the initial coordinates of the grid points on the lower surface of the top plate, and the initial coordinates of the grid points on the upper surface of the bottom plate. Within the displacement-limited area, the coordinates of the grid on the lower surface of the top plate are used as elements to generate a top plate grid list in ascending order of the corresponding grid ID in the numerical model. The elements in the top plate grid list are sorted, and the coordinates of the bottom plate grids that coincide with the horizontal plane of each top plate grid are sorted and used as elements to generate a bottom plate grid list. The coordinates of all grid points on the lower surface of the top plate within the displacement-limited area are sorted and used as elements to generate a top plate grid point list. The elements in the top plate grid point list are sorted, and the coordinates of the bottom plate grid points that coincide with the horizontal plane of each top plate grid point are sorted and used as elements to generate a bottom plate grid point list. Elements with the same ID in the top and bottom grid lists are treated as grid embedding pairs, and elements with the same ID in the top and bottom grid point lists are treated as grid point embedding pairs. The value of the ID is the number of list elements after the data is added. s4. Based on mesh embedding pairs and mesh point embedding pairs, a small deformation secondary simulation is performed on the initial numerical model. During the simulation, the numerical model information of the small deformation secondary simulation is extracted every 1-1 / 2 s steps, including the displacement vectors of the mesh points represented by the elements in the top and bottom plate mesh point lists. Each displacement vector is added as an element to the real-time top plate displacement list. According to the recording order of the top plate mesh points corresponding to each element in the real-time top plate displacement list, the order of the corresponding bottom plate mesh points is obtained according to the mesh point embedding pair relationship. The displacement parameters of each bottom plate mesh point are extracted according to the bottom plate mesh point order and added to the real-time bottom plate displacement list. The real-time surrounding rock deformation list is obtained by subtracting the corresponding element in the real-time bottom plate displacement list from each element in the real-time top plate displacement list. The positions of elements whose deformation components in the surrounding rock deformation list are greater than the dimension of this component in the excavated area are mapped to the real-time top and bottom plate displacement list. This yields the real-time displacement limit area in the process of calculating the unbalanced force to the target value using the numerical model of small deformation secondary simulation. Displacement correction is achieved by fixing the velocity of the large displacement side grid points in the top and bottom plates within the real-time displacement limit area to 0. The displacement correction in the real-time displacement limit area is removed using FLAC3D software commands, and an equivalent reaction force is applied to the grid points in the area where the displacement correction is removed. The reaction force of the grid points in the original real-time displacement limit area is extracted. Combining the embedding relationship between the top plate grid points and the bottom plate grid points, pressure is applied to the corresponding grid points in the top plate when extracting the bottom plate grid points, and pressure is applied to the corresponding grid points in the bottom plate when extracting the corresponding grid points in the top plate. After the pressure is applied, the displacement on the stress extraction side is fixed to 0. The mesh stress of the shallow layer of the self-excavated area on the displacement-limited side is extracted using FLAC3D software. The shallow layer is the first stratum outward from the top and bottom plates of the excavated area, and is set according to actual needs, using an approximate calculation formula: Transformation, combined with mesh embedding, corrects the mesh stress of the first layer of the bottom plate corresponding to the top plate or the bottom plate corresponding to the top plate from the excavated area outward, and generates displacement and stress correction functions between rock layers or surrounding rocks. Where, σ i ’ The stress in the shallow mesh on the side to be corrected is given by: σ0 is the stress in the shallow mesh on the side to be corrected as calculated in the simulation before correction; E2 is the Young's modulus of the side to be corrected; E1 is the Young's modulus of the shallow mesh on the displacement-constrained side; μ1 is the Poisson's ratio of the shallow mesh on the displacement-constrained side; and μ2 is the Poisson's ratio of the shallow mesh on the side to be corrected. j To simulate the stress in the shallow mesh under displacement constraints after a certain number of steps, σ j0 To simulate the stress in the shallow layer of the side-shallow mesh within a certain number of steps of pre-displacement, ε i The strain of the shallow mesh on the side to be corrected; s5. Repeat steps s2-s4 above until the grid point coordinates of the rock strata and surrounding rock are superimposed with displacement parameters and there is no mutual embedding.

2. The FLAC3D small deformation mode approximation method based on displacement and stress correction according to claim 1, characterized in that: The displacement-limited region is larger than the region where displacement parameters are interleaved after the coordinates of the small deformation simulation mesh are superimposed.

3. The FLAC3D small deformation mode approximation method based on displacement and stress correction according to claim 1, characterized in that: During the implementation of the method, displacement correction includes three processes: traversal, judgment and restriction, while stress correction includes three processes: reaction force extraction, pressure application and stress correction.

4. The FLAC3D small deformation mode approximation method based on displacement and stress correction according to claim 3, characterized in that: The displacement correction process iterates through the entire list of grid points and all elements of the grid list in the numerical model; the judgment process is determined by comparing the sum of the original coordinates and real-time displacement of the grid point with the actual position limit range; the constraint process is achieved by fixing the node velocity in one or more directions to 0 and deleting grid points in the list whose displacement has been constrained.

5. The FLAC3D small deformation mode approximation method based on displacement and stress correction according to claim 4, characterized in that, The method for obtaining the stress correction function is as follows: Reaction extraction involves removing the displacement constraint and applying an equivalent reaction force using the "zone gridpoint force-reaction" command, and then using a loop to traverse the grid points within the displacement-limited area to extract the reaction force of the nodes in the original constraint area. Pressure application involves calculating the approximate stress in the vicinity of each grid point by using the reaction force of the grid points on the displacement-fixed side and the area of ​​the adjacent region. Using the "zone face apply stress" command, the grid points on the non-displacement-fixed side are obtained through the grid point embedding relationship. Stress is applied to the grids that are in contact with the grid points on the non-displacement-fixed side. After the pressure application is completed, the displacement on the stress extraction side is fixed to 0. Stress correction involves extracting the stress from the shallow mesh of the top plate and applying an approximate calculation formula. The transformation, combined with mesh embedding, is a correction function for the mesh stress of the first layer of the bottom plate from the excavated region outwards, where σ i ’ Shallow stress in the base plate, σ0 is the uncorrected shallow stress in the base plate, E2 is the Young's modulus of the base plate, E1 is the Young's modulus of the top plate, μ1 is the Poisson's ratio of the top plate, μ2 is the Poisson's ratio of the base plate, σ j σ j To simulate the stress in the shallow top plate after a certain number of steps, σ j0 To simulate the stress in the shallow layer of the top plate before a certain number of steps, ε i This represents the shallow strain of the base plate.

6. The FLAC3D small deformation mode approximation method based on displacement and stress correction according to claim 5, characterized in that: After the numerical model of the small deformation secondary simulation is self-equilibrated, the "model step" command is used for solving. According to the parameters corresponding to the preset number of numerical simulation steps, the displacement correction function and stress correction function must be called after each line of solution command. When the numerical model is not balanced, the displacement correction function and stress correction function are called once after each preset running time.

7. A computer device, characterized in that, It includes a processor and a memory, the processor being electrically connected to the memory, the memory being used to store instructions and data, and the processor being used to execute the FLAC3D small deformation mode approximation simulation method based on displacement and stress correction as described in any one of claims 1-6.

Images