A Mesh Generation Method for Numerical Simulation of Curtain Grouting in Fractured Rock Mass

Abstract

The present invention relates to a grid meshing method for numerical simulation of curtain grouting in fractured rock masses. First, a geometric model of the fractured rock mass is constructed, which includes filled surface graphics, hollowed surface graphics, and line graphics. Then, according to the positions of the respective graphics, all intersecting graphics in the model are obtained. By traversing all the intersecting graphics, all the intersection points and the endpoints of the common edge segments are used as the two-dimensional fixed points of the model. Finally, global uniform point distribution is carried out within the model area to obtain the initial nodes. The initial nodes are screened according to the probability of the node positions to obtain the retained initial nodes, and the initial triangular mesh is generated from the retained initial nodes. The positions of the initial triangular mesh nodes are adjusted according to the spring principle. Considering one-dimensional constraints, based on all the nodes after position adjustment and the boundaries of all the graphics, the model is triangulated using the Delaunay triangulation method considering constraints, and the mesh quality is evaluated. This method takes into account one-dimensional and two-dimensional constraints during the grid meshing process, avoids the mutual intrusion of grid elements, and improves the grid quality.

Description

Technical Field

[0001] The present invention belongs to the technical field of numerical simulation of curtain grouting, and particularly relates to a mesh generation method for numerical simulation of curtain anti-seepage grouting in fractured rock masses. Background Art

[0002] Bedrock anti-seepage is an important part of projects such as the construction of water conservancy hubs, mine sewage treatment, and pumped storage power station construction. Curtain grouting is one of the most widely used bedrock anti-seepage measures, and numerical simulation is an important means to reveal the anti-seepage mechanism of grouting curtains, quantitatively evaluate the anti-seepage efficiency of curtains, and optimize the design of grouting schemes. Due to long-term geological effects, a large number of discontinuous structural planes appear in rock masses, such as faults, joints, and fractures. These discontinuous structural planes play a decisive role in the mechanical properties and seepage characteristics of rock masses, making the rock masses exhibit high anisotropy, inhomogeneity, and discontinuity.

[0003] In traditional bedrock grouting, the grouted medium is generally modeled using a continuous medium model. There is a large difference between the continuous medium model and the actual distribution of fractures. The grouting curtain formed by grouting is assumed to be continuous, complete, uniform, and completely impermeable. Therefore, it is impossible to simulate the flow of grout, and the interpenetration characteristics of grouted fractures and ungrouted fractures under the action of seepage filtration effect are ignored, resulting in inaccurate simulation of grouting effects. With the maturity of the technology for identifying drilling parameters while drilling and downhole imaging technology, the fractured rock mass has gradually changed from a "completely uncertain" system to a "partially known, partially unknown" system, making it possible to reconstruct the fracture network model of the grouted rock mass. The unified pipe network method proposed by Ma Guowei et al. in the article "Unified pipe network method for simulation of water flow in fractured porous rock" is increasingly used in the multi-field coupling numerical calculation of fractured pore media, enabling refined numerical simulation of the flow and diffusion of grout.

[0004] Constructing a grid model through numerical discretization is an important prerequisite for simulating slurry flow and forming a curtain to achieve refined simulation. Since the quality of the grid after model discretization is directly related to the computational complexity and accuracy, scholars have conducted extensive research on model grid meshing. The literature "A simple mesh generator in Matlab" established an unstructured grid meshing method through a spring model, but this method did not consider constrained grid meshing. For details, see Persson P O, Strang G. A simple mesh generator in Matlab[J], SIAM Review. 2004, 46(02): 329–345. The literature "Selective refinement: A new strategy for automatic node placement in graded triangular meshes" proposed a constrained arbitrary polygon Delaunay triangulation algorithm and an automatic triangular mesh node placement strategy. However, when dealing with models with a large number of fractures and complex spatial distributions, a large number of poor-quality grids will appear. For details, see Frey W H. Selective refinement: A new strategy for automatic node placement in graded triangular meshes[J], International Journal for Numerical Methods in Engineering. 1987, 24(11): 2183–2200.

[0005] In summary, the main problems of the existing technologies are as follows: (1) The fracture distribution in the model is quite different from the actual situation; (2) The grid meshing does not consider fracture constraints and ignores the potential impact of fracture uncertainty on the formation of leakage channels. Therefore, it is difficult to achieve refined simulation of curtain grouting in fractured rock masses, resulting in inaccurate prediction results. Summary of the Invention

[0006] Aiming at the deficiencies of the existing technologies, the technical problem to be solved by the present invention is to provide a grid meshing method for numerical simulation of curtain anti-seepage grouting in fractured rock masses.

[0007] The technical solution adopted by the present invention to solve the above technical problem is as follows:

[0008] A grid meshing method for numerical simulation of curtain anti-seepage grouting in fractured rock masses, characterized in that the method comprises the following steps:

[0009] Step 1: Construct a geometric model of fractured rock mass. In the model, the mountain body and water body are represented by filled areas, the grouting gallery and grouting boreholes are represented by hollowed areas, and geological faults and random joint fractures are represented by random line segments;

[0010] Step 2: Consider a closed filled area in the model as a filled surface graph, a closed hollowed area as a hollowed surface graph, and a random line segment as a line graph; According to the positions of each graph, obtain all the intersecting graphs in the model; Traverse all the intersecting graphs, and regard all the intersection points and the endpoints of the common-edge line segments as the two-dimensional fixed points of the model;

[0011] Step 3: Conduct triangular mesh generation for the model by the Delaunay triangulation method considering constraints;

[0012] Step 3.1: Uniformly arrange nodes on all line graphs with the minimum grid size as the node distance. The arranged nodes are the one-dimensional fixed points of the model; According to the connection relationship between the one-dimensional fixed points, obtain the one-dimensional constraints;

[0013] Step 3.2: Conduct global uniform point distribution within the model area with the minimum grid size as the spacing to obtain the initial nodes; Calculate the probability of the position of each node according to Equation (1);

[0014] Φ(p) = 1 / h(p).^2 (1)

[0015] h(p) = min(h 0 +h grad ·d(p),h 0 ) (2)

[0016]

[0017] where, Φ(p) represents the probability function of the node position, h(p) represents the node density function, d(p) represents the global distance function of the node, p represents the coordinate matrix of all nodes, h 0 represents the minimum grid size, h grad represents the grid size gradient, represents the boundary of the filled surface graph Ω i of, represents the boundary of the hollowed surface graph Π i of, represents the boundary of the line graph Γ i of, F j represents the jth two-dimensional fixed point;

[0018] For a certain node position, generate a random number between 0 and 1. If the random number is less than or equal to the probability of that node position, retain the initial node at that position; otherwise, delete it. Similarly, traverse all node positions, screen the initial nodes, and obtain the retained initial nodes;

[0019] Step 3.3: Generate an initial triangular mesh from the retained initial nodes according to the Delaunay triangulation method;

[0020] Step 3.4: Adjust the positions of the nodes of the initial triangular mesh according to the spring principle;

[0021] Step 3.5: Considering one-dimensional constraints, based on all the nodes after position adjustment and the boundaries of all the figures, use the Delaunay triangulation method considering constraints to perform triangular mesh generation on the model;

[0022] Step 4: Evaluate the quality of each triangular mesh according to Equation (6);

[0023]

[0024] Among them, q represents the quality score, r in represents the inradius of the triangular mesh, r out represents the circumradius of the triangular mesh, and a, b, and c respectively represent the three side lengths of the triangular mesh;

[0025] If the quality score of a certain triangular mesh is greater than or equal to the quality evaluation threshold, retain the triangular mesh; otherwise, delete all the edges of the triangular mesh, form a polygonal blank area in the area where the triangular mesh is located, and regenerate a triangular mesh in the polygonal blank area until the quality score of the regenerated triangular mesh is greater than or equal to the quality evaluation threshold.

[0026] Furthermore, in Step 3.4, calculate the repulsive force received by each node according to Equation (4);

[0027]

[0028] Among them, δ represents the spring coefficient, l kw represents the distance between the current node k and the neighbor node w, represents the ideal distance between the node k and the neighbor node w;

[0029] Adjust the positions of each node according to Equation (5):

[0030] k(t n+1 ) = k(t n ) + Δt·f(k) (5)

[0031] Among them, k(t n+1 )、k(tn ) represent the positions of node k at the (n + 1)-th and n-th iterations respectively, and t n+1 and t n represent the time steps at the (n + 1)-th and n-th iterations respectively, and Δt represents the time difference between the (n + 1)-th iteration and the n-th iteration;

[0032] Repeat the above operations to perform iterative solution on the positions of each node until the forces on each node are in balance, then stop the iteration; otherwise, return to step 3.3.

[0033] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0034] 1. In terms of model generation, the present invention introduces a complex fracture network, which not only considers large geological faults with limited quantity and determined occurrence, but also considers small joint fractures with huge quantity and relatively random spatial distribution. Compared with the traditional continuous medium model, the present invention considers the complex spatial structure of the grouted bedrock fracture network, can more realistically simulate the discontinuity and non-uniformity of the grouted rock mass, and provides a model basis for more accurate calculation of the anti-seepage efficiency of the curtain body. For fracture information, the present invention obtains data such as the occurrence, dip direction, and dip angle of fractures from geological exploration profiles, and generates a Stochastic Discrete Fracture Network based on these data.

[0035] 2. In terms of model discretization, the model grid meshing of the present invention is generated based on the CDT method, and parallel algorithms are used respectively when obtaining the one-dimensional fixed point and two-dimensional fixed point of the model, significantly improving the calculation efficiency and increasing the model processing capacity to thousands of fractures. At the same time, due to the grid meshing considering constraints strictly follows one-dimensional constraints and two-dimensional constraints, avoiding the mutual intrusion of grid cells, and the quality inspection after grid generation further corrects inferior grids, improving the grid quality. It is worth emphasizing that in the present invention, the one-dimensional constraint of the grid not only includes one-dimensional line graphics, but also includes the common edges of two-dimensional surface graphics, which makes the numerical calculation of the medium interface more convenient. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 is the overall flowchart of the present invention;

[0037] Figure 2 is the structural schematic diagram of the fractured rock mass geometric model of the present invention;

[0038] Figure 3 is the schematic diagram of the principle of the Delaunay triangulation method of the present invention;

[0039] Figure 4 is the schematic diagram of the grid meshing result of the present invention. DETAILED DESCRIPTION OF THE INVENTION

[0040] Specific embodiments are given below in conjunction with the accompanying drawings. The specific embodiments are only used to illustrate the technical solutions of the present invention in detail and are not used to limit the protection scope of this application.

[0041] The present invention is a mesh generation method for numerical simulation of curtain grouting in fractured rock masses (hereinafter referred to as the method, see Figures 1 to 4 ), which includes the following steps:

[0042] Step 1: Based on the grouting curtain structure design drawing and the geological exploration profile, use multi-physics field simulation software (Comsol) or programming software (Matlab, Python, etc.) to construct a geometric model of fractured rock masses; as Figure 1 shown, the fractured rock masses include mountains, water bodies, grouting galleries, grouting boreholes, geological faults, and random joint fractures. In the model, the mountains and water bodies are represented by filled areas, the grouting galleries and grouting boreholes are represented by hollowed-out areas, and the geological faults and random joint fractures are represented by discontinuous areas; among them, the filled areas are represented by a combination of polygons and ellipses, the discontinuous areas are represented by random line segments, and the hollowed-out areas are represented by any geometric figure other than random line segments.

[0043] Step 2: A closed filled area in the model is regarded as a filled surface figure, a closed hollowed-out area is regarded as a hollowed-out surface figure, and a random line segment is regarded as a line figure. Therefore, the model contains three types of figures: filled surface figures, hollowed-out surface figures, and line figures; according to the positional relationship of each figure in the geometric model of fractured rock masses, obtain the two-dimensional fixed points of the model.

[0044] All filled surface figures form a set Ω = {Ω 1 , Ω 2 , …, Ω i , …, Ω M}, all hollowed-out surface figures form a set Π = {Π 1 , Π 2 , …, Π i , …, Π M}, all line figures form a set Γ = {Γ 1 , Γ 2 , …, Γ i , …, Γ M}, i = 1, 2, …, M, where M represents the number of figures, and the number of each type of figure may be different;

[0045] For any filled surface figure Ω i , according to the positional relationship between the filled surface figure Ω i and each of the remaining figures (including the remaining filled surface figures, all hollowed-out surface figures, and line figures), obtain the relationship with the filled surface figure Ω iIntersecting figures; similarly, traverse all filled surface figures, hollowed surface figures, and line figures to obtain all intersecting figures; traverse all intersecting figures, and take all intersection points and the endpoints of the common edge segments as the two-dimensional fixed points of the model. All two-dimensional fixed points form a set F = {F 1 , F 2 , …, F j , …, F R}, where j = 1, 2, …, R, and R represents the number of two-dimensional fixed points;

[0046] Step 3: Perform global uniform point distribution within the model area to obtain initial nodes; screen the initial nodes according to the probability of the node positions to obtain the retained initial nodes; generate an initial triangular mesh from the retained nodes according to Delaunay triangulation; use an iterative solution method to adjust the positions of the initial triangular network nodes to obtain the nodes with adjusted positions; consider the one-dimensional constraints of the model, and use the Constrained Delaunay Triangulation (CDT) method to perform triangular mesh division on the model according to all the nodes with adjusted positions and the boundaries of all figures;

[0047] Step 3.1: For all line figures, evenly distribute nodes on the line figures with the minimum mesh size as the node distance. This node is the one-dimensional fixed point of the model, and at the same time, obtain one-dimensional constraints to represent the connection relationship between one-dimensional fixed points; the one-dimensional constraints are numerically 0 or 1, where 0 means two nodes are not connected, and 1 means two nodes are connected; for example, the two endpoints of a certain line figure are denoted as nodes u 1 and u 2 . Nodes u 1 and u 2 are arranged with node u 3 between them. Then the one-dimensional constraint formed by node u 1 and u 3 is 1, the one-dimensional constraint formed by node u 1 and u 2 is 0, and the one-dimensional constraint formed by node u 3 and u 2 is 1;

[0048] Step 3.2: Perform global uniform point distribution within the model area with the minimum mesh size as the spacing, and calculate the probability of each node position according to Equation (1);

[0049] Φ(p) = 1 / h(p).^2 (1)

[0050] h(p) = min(h 0 + h grad ·d(p), h 0 ) (2)

[0051]

[0052] Among them, Φ(p) represents the probability function of the node position, h(p) represents the node density function, d(p) represents the global distance function of the node, and p represents the coordinate matrix of all nodes. represents the filled surface graph Ω i of the boundary represents the hollowed surface graph Π i of the boundary represents the line graph Γ i of the boundary, h 0 represents the minimum grid size, h grad represents the grid size gradient;

[0053] The global distance function represents the minimum distance between any position in the model area and each graph. The density function represents the node distribution. From the place with a smaller distance function value to the place with a larger distance function value, the node distribution gradually becomes sparse, the node density decreases, and the grid size increases.

[0054] For a certain node position, a random number from 0 to 1 is generated. If the random number is less than or equal to the probability of the node position, the initial node at this node position is retained; otherwise, it is deleted. Similarly, traverse all node positions, and screen the initial nodes according to the probability of the node position to obtain the retained initial nodes. Assume that the probability of a certain node position is 0.8, and the probability that the initial node at this node position is retained is 80%. A random number from 0 to 1 is generated. If the random number is less than or equal to 0.8, the initial node at this node position is retained; if the random number is greater than 0.8, the initial node at this node position is deleted.

[0055] Step 3.3: According to the Delaunay Triangulation (DT) method, generate an initial triangular mesh from the retained initial nodes.

[0056] Step 3.4: Adjust the positions of the initial triangular mesh nodes according to the spring principle.

[0057] For the triangular mesh nodes, use the spring principle to adjust the positions of the nodes so that each node reaches a state of force balance, and then make the triangular mesh reach an ideal size. The spring principle is to regard the distance between adjacent nodes as a repulsive force, and the edges of the connected triangular meshes are considered as ideal springs, and adjust the node positions to make them reach a state of force balance.

[0058] Calculate the repulsive force received by each node according to Equation (4); for any node k, that is, calculate the repulsive force received by node k from all its neighbor nodes connected to it through the triangular mesh.

[0059]

[0060] Among them, δ represents the spring coefficient, and l kw represents the distance between the current node k and the neighbor node w, and represents the ideal distance between the node k and the neighbor node w;

[0061] Adjust the position of node k according to Equation (5) to make it reach the position of force balance;

[0062] k(t n+1 ) = k(t n ) + Δt·f(k) (5)

[0063] Among them, k(t n+1 ) and k(t n ) respectively represent the positions of node k in the (n + 1)-th and n-th iterations, t n+1 , t n respectively represent the time steps of the (n + 1)-th and n-th iterations, and Δt represents the time difference between the (n + 1)-th iteration and the n-th iteration;

[0064] Iteratively solve the positions of each node until each node reaches force balance, then stop the iteration and complete the position adjustment of the nodes; otherwise, return to step 3.3 to regenerate the initial triangular mesh;

[0065] Step 3.5: Considering one-dimensional constraints, according to all the nodes after position adjustment and the boundaries of all the figures, use the Delaunay triangulation method considering constraints to perform triangular mesh generation on the model;

[0066] Step 4: Evaluate the quality of each triangular mesh according to Equation (6) to obtain the quality score;

[0067]

[0068] Among them, q represents the quality score, r in represents the inradius of the triangular mesh, r out represents the circumradius of the triangular mesh, and a, b, c respectively represent the three side lengths of the triangular mesh;

[0069] If the quality score is greater than or equal to the quality assessment threshold, for example 0.5, it is considered that the quality of the triangular mesh meets the requirements, and the triangular mesh is retained; otherwise, all the edges of the triangular mesh are deleted. After deleting the triangular mesh, the remaining edges of the triangular meshes sharing edges with this triangular mesh together form a polygonal blank area. In this polygonal blank area, triangular meshes are regenerated, that is, points are evenly distributed in the polygonal blank area with the minimum mesh size as the node spacing, and triangular meshes are generated according to the Delaunay triangulation principle. The positions of the nodes are adjusted by the spring principle, and the quality of the regenerated triangular meshes is evaluated until the quality score of the triangular meshes is greater than or equal to the quality assessment threshold.

[0070] What is not described in the present invention is applicable to the prior art.

Claims

1. A grid meshing method for numerical simulation of curtain grouting for seepage prevention in fractured rock masses, characterized in that, the method comprises the following steps: Step 1: Construct a geometric model of the fractured rock mass. In the model, the mountain body and water body are represented by filled areas, the grouting gallery and grouting boreholes are represented by hollowed-out areas, and geological faults and random joint fractures are represented by random line segments; Step 2: Consider a closed filled area in the model as a filled surface figure, a closed hollowed-out area as a hollowed-out surface figure, and a random line segment as a line figure; according to the positions of each figure, obtain all the intersecting figures in the model; traverse all the intersecting figures, and regard all the intersection points and the endpoints of the common-edge line segments as the two-dimensional fixed points of the model; Step 3: Perform triangular grid meshing on the model by using the Delaunay triangulation method considering constraints; Step 3.1: Uniformly arrange nodes on all line figures with the minimum grid size as the node distance. The arranged nodes are the one-dimensional fixed points of the model; according to the connection relationship between the one-dimensional fixed points, obtain one-dimensional constraints; Step 3.2: Perform global uniform point distribution in the model area with the minimum grid size as the spacing to obtain initial nodes; calculate the probability of each node position according to formula (1); Φ(p) = 1 / h(p).^2 (1) h(p) = min(h 0 + h grad · d(p), h 0 ) (2) Among them, Φ(p) represents the probability function of the node position, h(p) represents the node density function, d(p) represents the global distance function of the node, p represents the coordinate matrix of all nodes, h 0 represents the minimum grid size, h grad represents the grid size gradient, represents the boundary of the filled surface graph Ω i ; represents the boundary of the hollowed surface graph Π i ; represents the boundary of the line graph Γ i ; F j represents the j-th two-dimensional fixed point; For a certain node position, generate a random number between 0 and 1. If the random number is less than or equal to the probability of this node position, retain the initial node at this node position; otherwise, delete it. Similarly, traverse all node positions and screen the initial nodes to obtain the retained initial nodes; Step 3.3: Generate an initial triangular grid from the retained initial nodes according to the Delaunay triangulation method; Step 3.4: Adjust the positions of the nodes of the initial triangular grid according to the spring principle; Step 3.5: Consider one-dimensional constraints, and perform triangular grid meshing on the model by using the Delaunay triangulation method considering constraints according to all the nodes after position adjustment and the boundaries of all figures; Step 4: Evaluate the quality of each triangular grid according to formula (6); Among them, q represents the quality score, and r in represents the inradius of the triangular mesh, and r out represents the circumradius of the triangular mesh, and a, b, and c respectively represent the three side lengths of the triangular mesh; If the quality score of a certain triangular grid is greater than or equal to the quality evaluation threshold, retain this triangular grid; otherwise, delete all the edges of this triangular grid, form a polygonal blank area in the area where this triangular grid is located, and regenerate a triangular grid in the polygonal blank area until the quality score of the regenerated triangular grid is greater than or equal to the quality evaluation threshold.

2. The grid meshing method for numerical simulation of curtain grouting for seepage prevention in fractured rock masses according to claim 1, characterized in that, in Step 3.4, calculate the repulsive force received by each node according to formula (4); where δ represents the spring coefficient, and l kw represents the distance between the current node k and its neighbor node w, and represents the ideal distance between node k and neighbor node w; Adjust the positions of each node according to formula (5): k(t n+1 ) = k(t n ) + Δt·f(k) (5) Among them, k(t n+1 ) and k(t n ) represent the positions of node k at the (n + 1)-th and n-th iterations respectively, t n+1 and t n represent the time steps of the (n + 1)-th and n-th iterations respectively, and Δt represents the time difference between the (n + 1)-th iteration and the n-th iteration; Repeat the above operations to perform iterative solution on the positions of each node until the forces on each node are balanced, then stop the iteration; otherwise, return to Step 3.3.

Images