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

Abstract

The invention discloses an FLAC3D small deformation mode approximate simulation method based on displacement and stress correction, and belongs to the field of numerical simulation. Performing preliminary simulation on a target working condition based on a small deformation mode of FLAC3D software; a displacement cloud picture is extracted, grids and node embedding pairs are judged, and a displacement limiting area is determined according to an area where displacement does not conform to reality; extracting pointers of nodes and grids in the region, writing the pointers into a list, processing the grids and node pointers which are mutually embedded at the boundary, and forming a grid and node embedding pair list; compiling a displacement and stress correction function in combination with an actual object displacement interval; secondary simulation is carried out on the target working condition, and displacement and stress in the small deformation mode simulation process are corrected; and repeating the steps until no mutual embedding exists after the coordinates of the grid points are superposed with the displacement parameters. According to the method, the modeling requirement can be reduced, the method is suitable for a complex numerical model, and the distribution conditions of stress, displacement, plastic zones and the like of actual engineering are well reflected.

Description

Technical Field

[0001] The invention belongs to the field of numerical simulation, and in particular relates to a FLAC3D small deformation mode approximate simulation method based on displacement and stress correction. Background Art

[0002] At present, scientific research often adopts the methods of theoretical analysis, indoor experiments, numerical simulation, and field practice. Among them, indoor experiments require professional equipment and have certain experimental costs, and field practice is large-scale and requires coordination and linkage from all aspects. Numerical simulation has the advantages of low research cost, easy result acquisition, multi-factor and multi-level analysis, and multi-condition analysis at engineering scale, and is favored by many scholars.

[0003] As a commonly used numerical simulation software in geotechnical engineering, FLAC3D is based on the finite difference method and has the advantages of accurate and reasonable calculation results, applicability to unstable processes, small memory usage, and fast calculation speed. The FLAC3D calculation process can be divided into two modes: large deformation and small deformation. In the large deformation mode, the grid node coordinates change to a balanced state in real time, and the calculation results are more accurate and more widely applicable, but the calculation cost is higher and the grid deformation may cause program errors; in the small deformation mode, the node coordinates remain unchanged, the calculation speed is fast, and the program is not prone to errors, but it is not suitable for simulating large deformations. Summary of the invention

[0004] In view of the shortcomings of the existing technology, a FLAC3D small deformation mode approximate simulation method based on displacement and stress correction is provided. The method has simple steps and is easy to use. By simulating the working conditions and correcting the displacement and stress in the small deformation mode simulation process, there is no mutual embedding area after the grid point coordinates are superimposed with the displacement parameters, which effectively avoids the errors caused by grid deformation. The simulation results of working conditions with large deformation are more realistic than general small deformation calculations.

[0005] To achieve the above technical objectives, the present invention discloses a FLAC3D small deformation mode approximate simulation method based on displacement and stress correction, which specifically includes the following steps:

[0006] s1. Establish the initial numerical model in which the grid points on the top and bottom plates of the mining face and the tunneling roadway overlap each other 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, the displacement changes of the rock strata and surrounding rocks in the mining face and tunneling are preliminarily simulated with small deformations. The total number of simulation steps is recorded as S;

[0007] s2. Extract the coordinates of the grid points before and after the displacement from the preliminary simulation results. The coordinates before the displacement are the initial coordinates, and the coordinates after the small deformation displacement are the final coordinates. The initial coordinates of the grid points generate displacement parameters due to mining. The final coordinates of the grid points on the lower surface of the roof and all the grid points on the bottom rock layer in the numerical model are calculated by the formula A'=A+D, where A' represents the final coordinates of the grid points, 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 the grid points on the lower surface of the roof and the bottom rock layer are drawn using the above formula. The final coordinates of the grid points on the bottom rock layer are above the lower surface of the roof, which is recorded as the bottom plate embedded area. Similarly, the final coordinates of the grid points on the upper surface of the bottom plate and all the grid points on the top rock layer are calculated by the formula A'=A+D, and all the grid points on the lower surface of the roof and the bottom rock layer are drawn according to the final coordinates.

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

[0009] s3. Based on the displacement-limited area obtained in step s2, the initial coordinates of the surface grids and grid points of the top and bottom plates of the mining face and the tunneling roadway in the displacement-limited area are extracted, including the initial coordinates of the grids on the lower surface of the top plate, the initial coordinates of the grids 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] The grid coordinates of the lower surface of the top plate in the displacement-limited area are used as elements to generate a top plate grid list according to the ID numbers of the corresponding grids in the numerical model from small to large. According to the order of the elements in the top plate grid list, the bottom plate grid coordinates that coincide with the horizontal planes 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 in the displacement-limited area are sorted and used as elements to generate a top plate grid point list. According to the order of the elements in the top plate grid point list, the coordinates of the bottom plate grid points that coincide with the horizontal planes of each top plate grid point are sorted and used as elements to generate a bottom plate grid point list;

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

[0012] s4. Based on the grid embedding pairs and grid point embedding pairs, the initial numerical model is simulated with small deformation twice. During the simulation process, the numerical model information of the small deformation secondary simulation is extracted every 1-1 / 2S steps, including the displacement vectors of the grid points represented by the elements in the top and bottom plate grid point lists. Each displacement vector is added as an element into the real-time top plate displacement list. According to the record order of the top plate grid points corresponding to each element in the real-time top plate displacement list, the order of the corresponding bottom plate grid points is obtained according to the relationship of the grid point embedding pairs. The displacement parameters of each bottom plate grid point are extracted according to the order of the bottom plate grid points, and the displacement parameters are added to the real-time bottom plate displacement list. The real-time surrounding rock deformation list is obtained by subtracting the elements in the real-time bottom plate displacement list from the elements in the real-time top plate displacement list. The element position records of the deformation components of each element in the surrounding rock deformation list are greater than the dimension of this component in the hollowing area in the direction, and are mapped to the real-time top and bottom plate displacement lists to obtain the real-time displacement limit area in the process of calculating the unbalanced force to the target value by the numerical model of the secondary simulation of small deformation. The displacement correction is achieved by fixing the speed of the grid points on the large displacement side of the top and bottom plates in the real-time displacement limit area to 0. The displacement correction of the real-time displacement limit area is released through the FLAC3D software command, and an equivalent reaction force is applied to the grid points in the displacement correction release area. The reaction force of the grid points in the original real-time displacement limit area is extracted. Combined with the embedding relationship between the top plate grid points and the bottom plate grid points, pressure is applied to the corresponding grid points of the top plate when extracting the bottom plate grid points, and pressure is applied to the corresponding grid points of the bottom plate when extracting the corresponding grid points of the top plate.

[0013] After the pressure is applied, the displacement on the stress extraction side is fixed to 0, and the grid stress of the shallow layer of the hollowing area on the displacement limit side is extracted by FLAC3D software. The shallow layer is the first layer of stratum from the top and bottom plates of the hollowing area to the outside. The thickness is set according to actual needs, and the approximate calculation formula is used: i ’ =σ0+(E2 / E1 / (1-μ2 2 )*(σ i -μ1*Σ(σ j -σ j0 ) / 2))+E2*dε i / dt conversion, combined with the grid embedding pair list, corrects the grid stress of the bottom plate corresponding to the top plate or the top plate corresponding to the bottom plate from the hollowed-out area to the first layer outside, and generates the displacement and stress correction function between rock layers or surrounding rocks;

[0014] Among them, σ i ’ The shallow grid stress on the side to be corrected, σ0 is the shallow grid stress on the side to be corrected calculated by simulation before correction, E2 is the Young's modulus on the side to be corrected, E1 is the Young's modulus of the shallow grid on the displacement-limited side, μ1 is the Poisson's ratio of the shallow grid on the displacement-limited side, μ2 is the Poisson's ratio of the shallow grid on the side to be corrected, σj , σ j To simulate the stress of the shallow mesh on the displacement-limited side after a certain number of steps, σ j0 To simulate the stress of the shallow layer of the shallow grid on the side of the displacement limit before a certain number of steps, ε i is the strain of the shallow grid on the side to be corrected;

[0015] s5. Repeat the above steps s2-s4 until the grid point coordinates of the rock layer and surrounding rock are not embedded in each other after the displacement parameters are superimposed.

[0016] Furthermore, the displacement-limited region is larger than the region where the small deformation simulation grid coordinates are embedded with each other after superimposing the displacement parameters.

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

[0018] Furthermore, the content traversed during the displacement correction process is the grid point list of the entire numerical model and all elements of the grid list; the judgment process is determined by comparing the sum of the original coordinates of the grid points and the real-time displacement with the actual position limit range; the restriction process is achieved by fixing one or several directions of the node velocity to 0, and deleting the grid points in the list with restricted displacement.

[0019] Furthermore, the method for obtaining the stress correction function is as follows: reaction force extraction includes releasing the displacement restriction and applying an equivalent reaction force through the "zone gridpoint force-reaction" command, traversing the grid points in the displacement-restricted area in combination with a loop statement, and extracting the reaction force of the nodes in the original restricted area; pressure application includes calculating the approximate stress of the adjacent area of ​​each grid point through the reaction force of the grid point on the fixed side of the displacement and the area of ​​the adjacent area, and using the "zone face apply stress" command to obtain the grid point on the non-displacement fixed side through the grid point embedding relationship, applying stress to the grid in contact with the grid point on the non-displacement fixed side, and re-fixing the displacement of the stress extraction side to 0 after the pressure is applied; stress correction includes extracting the shallow grid stress of the top plate, using the approximate calculation formula σ i ’ =σ0+(E2 / E1 / (1-μ2 2 )*(σ i -μ1*Σ(σ j -σ j0 ) / 2))+E2*dε i / dt conversion, combined with the grid embedding pair list to correct the grid stress correction function of the first layer of the bottom plate from the hollow area to the outside, where σ i ’The shallow stress of the bottom plate, σ0 is the shallow stress of the bottom plate calculated by simulation before correction, E2 is the Young's modulus of the bottom 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 bottom plate, σ j , σ j To simulate the stress of the shallow roof 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 is the shallow strain of the bottom plate.

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

[0021] A computer device comprises a processor and a memory, wherein the processor is electrically connected to the memory, the memory is used to store instructions and data, and the processor is used to execute a FLAC3D small deformation mode approximate simulation method based on displacement and stress correction.

[0022] Beneficial effects: The present invention provides an approximate simulation method of FLAC3D small deformation mode based on displacement and stress correction. First, based on the FLAC3D small deformation mode, a preliminary simulation is performed on the target working condition, the displacement cloud map is extracted, the mesh and node embedding pairs are judged, and the displacement limit area is determined according to the displacement that does not conform to the actual area; secondly, the pointers of the nodes and meshes inside the area are extracted, written into a list, and the grid and node pointers embedded in each other at the boundary are processed to form a list of mesh and node embedding pairs; then, combined with actual practice, displacement and stress correction functions are written, the target working condition is simulated twice, and the displacement and stress in the small deformation mode simulation process are corrected; finally, the above steps are repeated until there is basically no displacement that does not conform to the actual area. This method can avoid program errors caused by mesh deformation, and the simulation results for working conditions with larger deformations are more realistic than general small deformation calculations. Compared with the prior art, it has the following advantages:

[0023] 1) The present invention expands the application scope of the small deformation mode and the simulation results are more accurate.

[0024] 2) The present invention reduces the modeling requirements for simulating working conditions with large deformations.

[0025] 3) The present invention avoids program errors caused by mesh deformation during the simulation process. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] Figure 1 It is a flow chart of the FLAC3D small deformation mode approximate simulation method based on displacement and stress correction of the present invention.

[0027] Figure 2 It is a displacement cloud diagram in the z direction of the small deformation mode in the embodiment of the present invention.

[0028] Figure 3 It 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 z-direction displacement cloud diagram after displacement and stress correction in an embodiment of the present invention.

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

[0031] The embodiments of the present invention are further described below in conjunction with the accompanying drawings:

[0032] The present invention discloses a FLAC3D small deformation mode approximate simulation method based on displacement and stress correction:

[0033] s1. Preliminary simulation of the target working condition based on the small deformation mode of FLAC3D software;

[0034] s2. Based on the preliminary simulation results of the target working condition, the relationship between the top grid and the grid point pair is constructed: the displacement cloud map is extracted, the grid and node embedding pairs are determined, and the displacement limit area is determined according to the mutual embedding area after the grid point coordinates are superimposed with the displacement parameters, and the displacement limit area is larger than the mutual embedding area after the grid point coordinates are superimposed with the displacement parameters;

[0035] s3. Based on the displacement-defined region and mesh and node embedding pairs in step s2, extract the pointers of nodes and meshes inside the region, write them into a list, process the mesh and node pointers embedded in each other at the boundary, and form a list of mesh and node embedding pairs;

[0036] s4. Based on the pointer list in step s3 and in combination with the actual situation, write a displacement correction function including traversal, judgment and restriction and a stress correction function including reaction force extraction, pressure application and stress correction;

[0037] s5. Based on the displacement and stress correction function of step s4, after the target working condition is simulated twice and the model of the small deformation mode simulation process is self-balanced, the displacement and stress correction function is called after each line of the solution command. When the model is not balanced, the displacement and stress correction function is called once for correction after each calculation of 1-1 / 2S (S is the total number of preliminary simulation steps);

[0038] s6. Repeat the above steps s2-s6 until there is no mutually embedded area after the grid point coordinates are superimposed with the displacement parameters.

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

[0040] s1. A numerical model in which the grid points on the top and bottom plates of the mining face and the tunneling roadway overlap on the horizontal plane is established. Based on the FLAC3D small deformation mode, a preliminary simulation is performed on the displacement changes of the rock strata and surrounding rocks in the mining face and tunneling. The total number of simulation steps is recorded as S; the z-direction displacement cloud diagram of the small deformation mode is shown in 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, and calculate the final coordinates of the grid points on the lower surface of the roof and all the grid points on the bottom rock layer by 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 the grid points on the lower surface of the roof and the bottom rock layer according to the final coordinates, and the final coordinates of the grid points on the bottom rock layer are above the lower surface of the roof, which is recorded as the bottom plate embedded area. Calculate the final coordinates of the grid points on the upper surface of the bottom plate and all the grid points on the roof rock layer by the formula A'=A+D, and draw all the grid points on the lower surface of the roof and the bottom rock layer according to the final coordinates. When the top and bottom plate surfaces meet, if they are not in a balanced state, they will continue to move along their respective displacement directions to form an embedding. Therefore, the final coordinates of the top plate rock layer below the upper surface of the bottom plate are recorded as the top plate embedding area. The lower surface of the top plate and the upper surface of the bottom plate are drawn using the coordinates. The area where the lower surface of the top plate and the upper surface of the bottom plate intersect is the mutual embedding area between the top plate and the bottom plate. The coordinates of each point in the mutual embedding area between the top plate and the bottom plate are obtained, and the three-axis data of the coordinates of each point in the mutual embedding area between the top plate and the bottom plate are sorted by size to obtain the maximum and minimum values ​​of the three-axis coordinates. The area between the maximum values ​​of the three-axis coordinates is determined as the displacement limit area; the displacement limit area is greater than the mutual embedding area after the small deformation simulation grid coordinates are superimposed on the displacement parameters;

[0042] s3. Based on the displacement limiting area obtained in step s2, extract the initial coordinates of the surface grids and grid points of the top and bottom plates of the mining face and the tunneling tunnel in the displacement limiting area, add all the grid coordinates of the lower surface of the top plate in the displacement limiting area into the top plate grid list, set the order of adding the grid coordinates of the lower surface of the top plate to the top plate grid list as the element order of the top plate grid point list, add the bottom plate grid coordinates that coincide with the horizontal planes of each top plate grid into the bottom plate grid list according to the element order in the top plate grid list, add all the grid point coordinates of the lower surface of the top plate in the area into the top plate grid point list, add the bottom plate grid point coordinates that coincide with the horizontal planes of each top plate grid point into the bottom plate grid point list according to the element order in 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, and 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 grid and grid point embedding pairs in step s3, perform a secondary step-by-step simulation of small deformation on the original model in step s1. During the secondary step-by-step simulation, extract the displacement parameters of the grid points represented by the elements in the top plate grid point list in the model of the secondary simulation of small deformation every 1-1 / 2S steps, and construct a real-time top plate displacement list for each displacement parameter. The displacement parameters in the real-time top plate displacement list correspond to the order of the top plate grid points: the displacement parameters in the real-time top plate displacement list exist in the order a according to the id number (addition time), and a group of displacement parameters (x, y, z) is generated from a top plate grid point b, so the top plate grid points also have an order a1, one top plate grid point corresponds to one bottom plate grid point b*, and when processing the bottom plate grid points, it is based on the order a of b. 1. Process b*; obtain the order of corresponding bottom plate grid points according to the embedding relationship of grid points, extract the displacement parameters of each bottom plate grid point according to the order of bottom plate grid points, add the displacement parameters to the real-time bottom plate displacement list, subtract the elements in the corresponding position in the real-time bottom plate displacement list from each element in the real-time top plate displacement list to obtain the real-time surrounding rock deformation list, map the element position record of the deformation component of each element in the real-time surrounding rock deformation list that is greater than the dimension element position record of this component direction of the hollowing area to the real-time top and bottom plate displacement list, and obtain the real-time displacement limit area in the process of calculating the unbalanced force to the target value in the model of the small deformation secondary simulation, and realize the displacement correction by fixing the speed of the large displacement side of the top and bottom plates in the real-time displacement limit area to 0; after the numerical model of the secondary simulation is self-balanced, use the "model step" command to solve, according to the parameters corresponding to the preset number of simulation steps, after each line of the solution command, the displacement and stress correction function must be called, and when the model is not balanced, the displacement and stress correction function is called once after each preset running time;

[0044] The displacement restriction of the real-time displacement limited area is released and the equivalent reaction force is applied through the FLAC3D software command, the reaction force of the original real-time displacement limited area node is extracted, and the embedding relationship between the top plate grid point and the bottom plate grid point is combined. When the bottom plate grid point is extracted, pressure is applied to the corresponding grid point of the top plate, and when the corresponding grid point of the top plate is extracted, pressure is applied to the corresponding grid point of the bottom plate. After the pressure is applied, the displacement of the stress extraction side is fixed to 0, and the shallow grid stress of the displacement limited side is extracted. The shallow layer is the first layer of stratum from the top and bottom plates of the hollowed area to the outside, and the thickness is set according to actual needs; the approximate calculation formula is used: σ i ’ =σ0+(E2 / E1 / (1-μ2 2 )*(σ i -μ1*Σ(σ j -σ j0 ) / 2))+E2*dε i / dt conversion, combined with the grid embedding pair list to correct the stress of the bottom plate corresponding to the top plate or the shallow grid of the top plate corresponding to the bottom plate, to generate the displacement and stress correction function between the rock layers or surrounding rocks;

[0045] In FLAC3D software, the displacement correction function includes three parts: traversal, judgment and restriction, and the stress correction function includes three parts: reaction force extraction, pressure application and stress correction; the displacement correction function traverses all elements of the node and grid pointer list; the judgment part makes a judgment by comparing the sum of the original coordinates of the node and the real-time displacement with the actual position limit range; the restriction part is achieved by fixing one or several directions of the node speed to 0, and deleting the node pointers with restricted displacement in the list.

[0046] Method for obtaining stress correction function: reaction force extraction includes releasing displacement restriction and applying equivalent reaction force through the "zone gridpoint force-reaction" command, traversing the node pair list in combination with loop statements, and extracting the reaction force of the nodes in the original restricted area; pressure application includes calculating the approximate stress of the adjacent area of ​​each node through the node reaction force and the area of ​​the adjacent area, and applying stress to the adjacent area on the contact side of the original node\grid by using the "zone face apply stress" command in combination with the embedded pair list. After the pressure is applied, the displacement of the stress extraction side is fixed to 0; stress correction includes extracting the shallow grid stress of the top plate, using the approximate calculation formula σ i ’ =σ0+(E2 / E1 / (1-μ2 2 )*(σ i -μ1*Σ(σ j -σ j0 ) / 2))+E2*dε i / dt conversion, combined with the grid embedding pair list to correct the correction function of the shallow grid stress of the bottom plate, where σi ’ The shallow stress of the bottom plate, σ0 is the shallow stress of the bottom plate calculated by simulation before correction, E2 is the Young's modulus of the bottom 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 bottom plate, σ j , σ j To simulate the stress of the shallow roof 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 is the shallow strain of the bottom plate.

[0047] Among them, σ i ’ The shallow grid stress on the side to be corrected, σ0 is the shallow grid stress on the side to be corrected calculated by simulation before correction, E2 is the Young's modulus on the side to be corrected, E1 is the Young's modulus of the shallow grid on the displacement-limited side, μ1 is the Poisson's ratio of the shallow grid on the displacement-limited side, μ2 is the Poisson's ratio of the shallow grid on the side to be corrected, σ j , σ j To simulate the stress of the shallow mesh on the displacement-limited side after a certain number of steps, σ j0 To simulate the stress of the shallow layer of the shallow grid on the side of the displacement limit before a certain number of steps, ε i is the strain of the shallow grid on the side to be corrected;

[0048] s5. Repeat the above steps s2-s4 until the grid point coordinates of the rock layer and surrounding rock are superimposed with displacement parameters without mutual embedding; effectively reflect the distribution of stress, displacement, plastic zone, etc. of the actual project, such as Figure 3 , Figure 4 and Figure 5 shown.

Claims

1. A FLAC3D small deformation mode approximate simulation method based on displacement and stress correction, characterized in that: The specific steps include: s1. Establish the initial numerical model in which the grid points on the top and bottom plates of the mining face and the tunneling roadway overlap each other 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, the displacement changes of the rock strata and surrounding rocks in the mining face and tunneling are preliminarily simulated with small deformations. The total number of simulation steps is recorded as S; s2. Extract the coordinates of the grid points before and after displacement from the preliminary simulation results. The initial coordinates are before displacement, and the final coordinates are after small deformation displacement. The initial coordinates of the grid points generate displacement parameters due to mining. The final coordinates of the grid points on the lower surface of the roof and all the grid points on the floor rock layer in the numerical model are calculated by the formula A'=A+D, where A' represents the final coordinates of the grid points, 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 layer are plotted using the above formula. The final coordinates of the grid points on the bottom rock layer are above the lower surface of the roof, which is recorded as the bottom plate embedded area. Similarly, the final coordinates of the upper surface grid points of the bottom plate and all grid points of the top plate rock layer are calculated by the formula A'=A+D, and all grid points of the lower surface of the top plate and the bottom plate rock layer are drawn according to the final coordinates; If the lower surface of the roof and the upper surface of the bottom plate are not in a balanced state when they meet in the small deformation simulation, the two will continue to move along their respective displacement directions to form an embedding. Therefore, the final coordinates of all grid points of the roof rock layer below the upper surface of the bottom plate are recorded as the roof embedding area; the lower surface of the roof and the upper surface of the bottom plate are drawn using the final coordinates of all grid points of the roof and bottom plate rock layers. The area where the lower surface of the roof and the upper surface of the bottom plate intersect is the top and bottom plate mutual embedding area. The coordinates of each grid point in the top and bottom plate mutual embedding area are obtained, and the three-axis data of the coordinates of each grid point in the top and bottom plate mutual embedding area are sorted by size to obtain the maximum and minimum values ​​of the three-axis coordinates, and the area between the maximum values ​​of the three-axis coordinates is determined as the displacement limit area; s3. Based on the displacement-limited area obtained in step s2, the initial coordinates of the surface grids and grid points of the top and bottom plates of the mining face and the tunneling roadway in the displacement-limited area are extracted, including the initial coordinates of the grids on the lower surface of the top plate, the initial coordinates of the grids 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; The grid coordinates of the lower surface of the top plate in the displacement-limited area are used as elements to generate a top plate grid list according to the ID numbers of the corresponding grids in the numerical model from small to large. According to the order of the elements in the top plate grid list, the bottom plate grid coordinates that coincide with the horizontal planes 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 in the displacement-limited area are sorted and used as elements to generate a top plate grid point list. According to the order of the elements in the top plate grid point list, the coordinates of the bottom plate grid points that coincide with the horizontal planes of each top plate grid point are sorted and used as elements to generate a bottom plate grid point list; The elements with the same id number in the top plate grid list and the bottom plate grid list are taken as grid embedding pairs, and the elements with the same id number in the top plate grid point list and the bottom plate grid point list are taken as grid point embedding pairs. The value of the id number is the number of list elements after the data is added. s4. Based on the grid embedding pairs and grid point embedding pairs, the initial numerical model is simulated with small deformation twice. During the simulation process, the numerical model information of the small deformation secondary simulation is extracted every 1-1 / 2S steps, including the displacement vectors of the grid points represented by the elements in the top and bottom plate grid point lists. Each displacement vector is added as an element into the real-time top plate displacement list. According to the record order of the top plate grid points corresponding to each element in the real-time top plate displacement list, the order of the corresponding bottom plate grid points is obtained according to the relationship of the grid point embedding pairs. The displacement parameters of each bottom plate grid point are extracted according to the order of the bottom plate grid points, and the displacement parameters are added to the real-time bottom plate displacement list. The real-time surrounding rock deformation list is obtained by subtracting the elements in the real-time bottom plate displacement list from the elements in the real-time top plate displacement list. The element position records of the deformation components of each element in the surrounding rock deformation list are greater than the dimension of this component in the hollowing area in the direction, and are mapped to the real-time top and bottom plate displacement lists to obtain the real-time displacement limit area in the process of calculating the unbalanced force to the target value by the numerical model of the secondary simulation of small deformation. The displacement correction is achieved by fixing the speed of the grid points on the large displacement side of the top and bottom plates in the real-time displacement limit area to 0. The displacement correction of the real-time displacement limit area is released through the FLAC3D software command, and an equivalent reaction force is applied to the grid points in the displacement correction release area. The reaction force of the grid points in the original real-time displacement limit area is extracted. Combined with the embedding relationship between the top plate grid points and the bottom plate grid points, pressure is applied to the corresponding grid points of the top plate when extracting the bottom plate grid points, and pressure is applied to the corresponding grid points of the bottom plate when extracting the corresponding grid points of the top plate. After the pressure is applied, the displacement on the stress extraction side is fixed to 0, and the grid stress of the shallow layer of the hollowing area on the displacement limit side is extracted by FLAC3D software. The shallow layer is the first layer of stratum from the top and bottom plates of the hollowing area to the outside. It is set according to actual needs and the approximate calculation formula is used: i ’ =σ0+(E2 / E1 / (1-μ2 2 )*(σ i -μ1*Σ(σ j -σ j0 ) / 2))+E2*dε i / dt conversion, combined with the grid embedding pair list, corrects the grid stress of the bottom plate corresponding to the top plate or the top plate corresponding to the bottom plate from the hollowed-out area to the first layer outside, and generates the displacement and stress correction function between rock layers or surrounding rocks; Among them, σ i ’ The shallow grid stress on the side to be corrected, σ0 is the shallow grid stress on the side to be corrected calculated by simulation before correction, E2 is the Young's modulus on the side to be corrected, E1 is the Young's modulus of the shallow grid on the displacement-limited side, μ1 is the Poisson's ratio of the shallow grid on the displacement-limited side, μ2 is the Poisson's ratio of the shallow grid on the side to be corrected, σ j , σ j To simulate the stress of the shallow mesh on the displacement-limited side after a certain number of steps, σ j0 To simulate the stress of the shallow layer of the shallow grid on the side of the displacement limit before a certain number of steps, ε i is the strain of the shallow grid on the side to be corrected; s5. Repeat the above steps s2-s4 until the grid point coordinates of the rock layer and surrounding rock are not embedded in each other after the displacement parameters are superimposed.

2. The FLAC3D small deformation mode approximate simulation method based on displacement and stress correction according to claim 1 is characterized in that: The displacement-limited area is larger than the area where the small deformation simulation grid coordinates are embedded with each other after superimposing the displacement parameters.

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

4. The FLAC3D small deformation mode approximate simulation method based on displacement and stress correction according to claim 3 is characterized in that: The content traversed during the displacement correction process is the grid point list of the entire numerical model and all elements of the grid list; the judgment process is determined by comparing the sum of the original coordinates of the grid point and the real-time displacement with the actual position limit range; the restriction process is achieved by fixing one or several directions of the node velocity to 0, and deleting the grid points in the list with restricted displacement.

5. The FLAC3D small deformation mode approximate simulation method based on displacement and stress correction according to claim 4 is characterized in that: Method for obtaining stress correction function: reaction force extraction includes releasing displacement restriction and applying equivalent reaction force through "zone gridpoint force-reaction" command, traversing grid points in displacement restricted area in combination with loop statement, and extracting reaction force of nodes in original restricted area; pressure application includes calculating approximate stress of adjacent area of ​​each grid point through reaction force of grid points on displacement fixed side and area of ​​adjacent area, and using "zone face apply stress" command to obtain grid points on non-displacement fixed side through grid point embedding relationship, and applying stress to grids in contact with grid points on non-displacement fixed side. After pressure application, the displacement of stress extraction side is fixed again to 0; Stress correction involves extracting the shallow grid stress of the roof and using the approximate calculation formula σ i ’ =σ0+(E2 / E1 / (1-μ2 2 )*(σ i -μ1*Σ(σ j -σ j0 ) / 2))+E2*dε i / dt conversion, combined with the grid embedding pair list to correct the grid stress correction function of the bottom plate from the hollow area to the first layer outside, where σ i ’ The shallow stress of the bottom plate, σ0 is the shallow stress of the bottom plate calculated by simulation before correction, E2 is the Young's modulus of the bottom 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 bottom plate, σ j , σ j To simulate the stress of the shallow roof 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 is the shallow strain of the bottom plate.

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

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

Images