An embedded finite element mesh generation method for three-dimensional fracture-matrix model

Abstract

The application discloses a kind of embedded finite element grid generation methods of three-dimensional fissure-matrix grid, comprising the following steps: 1) given matrix three-dimensional finite element grid, according to fissure space position and direction, find all intersection points of fissure surface and three-dimensional grid element;2) all intersection points are numbered in anticlockwise order;3) according to numbering order, form fissure element grid, record the three-dimensional matrix grid element number intersected with fissure element;4) check the quality of embedded fissure finite element grid, the application provides an accurate and efficient three-dimensional embedded fissure-matrix grid division method, solves the fissure grid modeling problem of embedded three-dimensional matrix model, and can be used to simulate the seepage and mechanical process of fissure and rock-soil coupling in rock-soil body containing fissure.

Description

Technical Field

[0001] This invention relates to the field of finite element modeling technology for fracture-pore coupled seepage simulation in the fields of hydraulic engineering and geotechnical engineering, specifically to an embedded finite element mesh generation method for a three-dimensional fracture-matrix model. Background Technology

[0002] One of the most prominent problems in water conservancy engineering construction is slope instability. Seepage is a key factor affecting slope stability. Before construction, appropriate seepage models, parameters, and boundary conditions are typically selected, and finite element analysis (FEM) is used for numerical simulation of seepage. This greatly helps in ensuring the safety of water conservancy projects, saving costs, and shortening time. When using finite element analysis for numerical simulation, matrix seepage under the influence of fractures must be considered. The quality of the three-dimensional finite element mesh determines the accuracy and speed of numerical calculations. However, the three-dimensional finite element mesh generation for fracture-embedded matrix models under complex geological conditions is extremely difficult. Using general commercial software to process complex models not only consumes a lot of manpower and time but also cannot effectively control the number and quality of meshes, resulting in a large computational load and poor simulation results. Summary of the Invention

[0003] The purpose of this invention is to overcome the above-mentioned shortcomings and provide an embedded finite element mesh generation method for a three-dimensional fracture-matrix model, which can be used for finite element analysis and calculation to achieve high-precision and high-efficiency fracture mesh generation and improve simulation calculation results.

[0004] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: an embedded finite element mesh generation method for a three-dimensional fracture-matrix model, which includes the following steps: 1) Given a matrix 3D finite element mesh, determine the intersection between the fracture surface and the generated matrix 3D finite element mesh based on the spatial location and orientation of the fracture, and find all the intersection points between them; 2) Determine the coordinates of the centroid of the polygon formed by all intersection points, and number all intersection points in counterclockwise order; 3) Form a fracture element mesh according to the numbering order, and record the three-dimensional matrix mesh element numbers that intersect with the fracture elements; 4) Check the quality of the embedded crack finite element mesh.

[0005] Further, step 1) includes: using the coordinate rotation method to find the intersection point between the fracture surface and the three-dimensional finite element mesh of the matrix, the specific steps of which are as follows: a. Let the plane containing any two intersecting edges of the fracture plane be the XOY plane, and let the endpoints A of the two edges be... B C Point A is the origin of the coordinate system, A→B is the positive X-axis, and the direction perpendicular to the XOY plane is the Z-axis. A new three-dimensional coordinate system is established. b. The coordinates of the finite element mesh nodes are changed from the original coordinate system. To the new coordinate system The transformation formula is:

[0006] In the formula, the rotation matrix R can be expressed as:

[0007] Matrix elements , , Let them be three pairwise orthogonal unit vectors, specifically represented as follows: set up It is the normal vector of the crack plane along the positive Z-axis, starting from point A. Recorded as , Recorded as Therefore , Then we have:

[0008] because , can be , Represented as:

[0009] Will Normalization yields:

[0010]

[0011] Finally, find the direction vector of the Y-axis:

[0012] Then Normalization yields:

[0013] If let middle, for , for , for ;but

[0014] c. Determine whether any of the 12 edges of the matrix element mesh intersect the fracture surface. In the new coordinate system, iterate through the 12 edges of the matrix mesh element. If the first edge intersects the fracture surface... The coordinates of the two endpoints of the edge are respectively and ,and Then the line segment must intersect the crack plane. Let the coordinates of the intersection point be... ,in , The calculation formula is:

[0015]

[0016] Further, step 2) includes: if there are There are several intersection points, and their coordinates in the new coordinate system are as follows: , , ..., The coordinates of the center point of the enclosed planar figure are The calculation formula is:

[0017]

[0018] The distance from any intersection point to the center point is The x-coordinate of the intersection point relative to the center point is Then the angle between the vector from the center point to the intersection point and the positive x-axis is:

[0019] If the ordinate of the intersection point is less than the ordinate of the center point Then the angle between it and the positive x-axis should be . Number the intersection points of all vectors with the positive x-axis in ascending order from 1, 2, ..., n, which is the counterclockwise order.

[0020] Further, step 3) includes: generating a fracture element mesh by sequentially generating the generated intersection point numbers, with a maximum of 8 intersection points. If the number of intersection points is less than 2, an element mesh cannot be formed; considering 6 cases where 3-8 intersection points are used to form a polygon to divide the fracture element mesh, and recording the 3D matrix mesh element number that intersects with the fracture element. If there are 3 intersection points, then the intersection points within the matrix grid cells will be combined into 1 fracture element grid, with the element numbering order being: 1→2→3→1; If the number of intersections is 4, then the intersections within the matrix grid cells will be combined into one fracture cell grid, with the cell numbering order being: 1→2→3→4; If the number of intersections is 5, then the intersections within the matrix mesh element will be combined into 2 fracture element meshes, with the element numbering order as follows: 1→2→3→4, 4→5→1→1; If the number of intersections is 6, then the intersections within the matrix grid cells will be combined into 2 fracture element grids, with the element numbering order as follows: 1→2→3→4, 4→5→6→1; If the number of intersections is 7, then the intersections within the matrix grid cells will be combined into 3 fracture element grids, with the element numbering order as follows: 1→2→3→4, 4→5→6→1, 6→7→1→1; If the number of intersections is 8, then the intersections within the matrix grid cells will be combined into 3 fracture element grids, with the element numbering order as follows: 1→2→3→4, 4→5→6→1, 6→7→8→1.

[0021] Furthermore, step 3) further includes: calculating the coordinate values ​​of the fracture mesh nodes in the original coordinate system. The calculation formula is: .

[0022] Furthermore, step 4) includes: checking the quality of the finite element mesh. The crack element mesh generated by the program may have individual elements with a Jacobian determinant of less than 0. The mesh quality is improved by manually adjusting the mesh shape.

[0023] The beneficial effects of this invention are: 1. This invention can be used to simulate seepage stress coupled finite element calculations in soil and rock masses containing discontinuous structural surfaces such as cracks and contact surfaces, and achieve high-efficiency and high-precision modeling of three-dimensional matrix meshes and embedded crack meshes, providing finite element models for numerical simulation calculations.

[0024] 2. Traditional methods typically model soil and rock media and fractures in a unified manner, making it difficult to achieve a transition between matrix and fracture meshes on a finite element basis. This results in a large finite element mesh size, impacting numerical computation efficiency. Embedded meshes treat fracture and matrix meshes as relatively independent systems, modeling them separately and achieving a mapping between them in the spatial domain. This eliminates the need to consider the transition between fracture and matrix meshes, significantly reducing the mesh size and improving numerical computation efficiency.

[0025] 3. Modeling the mapping relationship between the fracture mesh and the matrix 3D mesh is quite difficult. This invention proposes a calculation method based on the intersection of surfaces and line segments rotated in 3D space. This method obtains the intersection mode of the 3D fracture surface and the matrix 3D mesh element and establishes the mapping relationship between the fracture mesh and the matrix 3D mesh in the spatial domain. This method has high computational efficiency and reliable results.

[0026] 4. This invention cleverly utilizes the centroid of the intersection point between the fracture surface and the element, and determines the position of each intersection point by the angle between the line connecting each intersection point and the centroid and the local coordinate X-axis. It also realizes the sorting and numbering of fracture mesh nodes according to certain rules. This method is simple in principle and highly operable in three-dimensional space. Attached Figure Description

[0027] Figure 1 This is a flowchart illustrating the implementation process of the finite element mesh generation method for cracks in this invention. Figure 2 It is a new three-dimensional coordinate system diagram; Figure 3 It is a diagram showing the intersection of the matrix element mesh and the fracture plane; Figure 4 This is a diagram showing the intersection point numbers; Figure 5 This is a schematic diagram of the fracture element mesh generation; Figure 6 It is a three-dimensional mesh diagram of a given matrix; Figure 7 It is a planar mesh diagram of the fracture; Figure 8 It is a matrix three-dimensional finite element mesh generated by traditional methods; Figure 9 It is a three-dimensional finite element mesh of cracks generated by traditional methods. Detailed Implementation

[0028] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0029] like Figure 1 As shown, an embedded finite element mesh generation method for a three-dimensional fracture-matrix model includes the following steps: 1) Given a matrix 3D finite element mesh, determine the intersection between the fracture surface and the generated matrix 3D finite element mesh based on the spatial location and orientation of the fracture, and find all the intersection points between them; 2) Determine the coordinates of the centroid of the polygon formed by all intersection points, and number all intersection points in counterclockwise order; 3) Form a fracture element mesh according to the numbering order, and record the three-dimensional matrix mesh element numbers that intersect with the fracture elements; 4) Check the quality of the embedded crack finite element mesh.

[0030] Further, step 1) includes: using the coordinate rotation method to find the intersection point between the fracture surface and the three-dimensional finite element mesh of the matrix, the specific steps of which are as follows: a. Let the plane containing any two intersecting edges of the fracture plane be the XOY plane, and let the endpoints A of the two edges be... B C Point A is the origin of the coordinate system, A→B is the positive X-axis, and the direction perpendicular to the XOY plane is the Z-axis. A new three-dimensional coordinate system is established. b. The coordinates of the finite element mesh nodes are changed from the original coordinate system. To the new coordinate system The transformation formula is:

[0031] In the formula, the rotation matrix R can be expressed as:

[0032] Matrix elements , , Let them be three pairwise orthogonal unit vectors, specifically represented as follows: set up It is the normal vector of the fracture plane along the positive Z-axis, starting from point A. Recorded as , Recorded as Therefore , Then we have:

[0033] because , can be , Represented as:

[0034] Will Normalization yields:

[0035]

[0036] Finally, find the direction vector of the Y-axis:

[0037] Then Normalization yields:

[0038] If let middle, for , for , for ;but

[0039] c. Determine whether any of the 12 edges of the matrix element mesh intersect the fracture surface. In the new coordinate system, iterate through the 12 edges of the matrix mesh element. If the first edge intersects the fracture surface... The coordinates of the two endpoints of the edge are respectively and ,and Then the line segment must intersect the crack plane. Let the coordinates of the intersection point be... ,in , The calculation formula is:

[0040]

[0041] Further, step 2) includes: if there are There are several intersection points, and the coordinates of these intersection points in the new coordinate system are as follows: , , ..., The coordinates of the center point of the enclosed planar figure are The calculation formula is:

[0042]

[0043] The distance from any intersection point to the center point is The x-coordinate of the intersection point relative to the center point is Then the angle between the vector from the center point to the intersection point and the positive x-axis is:

[0044] If the ordinate of the intersection point is less than the ordinate of the center point Then the angle between it and the positive x-axis should be . Number the intersection points of all vectors with the positive x-axis in ascending order from 1, 2, ..., n, which is the counterclockwise order.

[0045] Further, step 3) includes: generating a fracture element mesh by sequentially generating the generated intersection point numbers, with a maximum of 8 intersection points. If the number of intersection points is less than 2, an element mesh cannot be formed; considering 6 cases where 3-8 intersection points are used to form a polygon to divide the fracture element mesh, and recording the 3D matrix mesh element number that intersects with the fracture element. If there are 3 intersection points, then the intersection points within the matrix grid cells will be combined into 1 fracture element grid, with the element numbering order being: 1→2→3→1; If the number of intersections is 4, then the intersections within the matrix grid cells will be combined into one fracture cell grid, with the cell numbering order being: 1→2→3→4; If the number of intersections is 5, then the intersections within the matrix mesh element will be combined into 2 fracture element meshes, with the element numbering order as follows: 1→2→3→4, 4→5→1→1; If the number of intersections is 6, then the intersections within the matrix grid cells will be combined into 2 fracture element grids, with the element numbering order as follows: 1→2→3→4, 4→5→6→1; If the number of intersections is 7, then the intersections within the matrix grid cells will be combined into 3 fracture element grids, with the element numbering order as follows: 1→2→3→4, 4→5→6→1, 6→7→1→1; If the number of intersections is 8, then the intersections within the matrix grid cells will be combined into 3 fracture element grids, with the element numbering order as follows: 1→2→3→4, 4→5→6→1, 6→7→8→1.

[0046] Furthermore, step 3) further includes: calculating the coordinate values ​​of the fracture mesh nodes in the original coordinate system. The calculation formula is: .

[0047] Furthermore, step 4) includes: checking the quality of the finite element mesh. The crack element mesh generated by the program may have individual elements with a Jacobian determinant of less than 0. The mesh quality is improved by manually adjusting the mesh shape.

[0048] Example 1: An embedded finite element mesh generation method for a three-dimensional fracture-matrix model includes the following steps: 1) The given matrix 3D model generated a total of 83,386 nodes and 79,935 elements, as detailed in the appendix. Figure 6 The coordinate system used to establish the 3D matrix model is rotated to form a new 3D coordinate system with the crack plane as the XOY plane, vertex A of the plane as the origin, and edge AB as the positive X-axis, as shown in the attached figure. Figure 2 As shown. Then the original coordinates Change to new coordinates The transformation is as follows:

[0049] Next, we need to determine the intersection of the 12 edges of the matrix element mesh with the fracture surface. First, we need to determine whether there are any intersection points. (See attached...) Figure 3 The mid-plane is a local fracture plane, and line segment ab is an edge of the matrix element hexahedron. If the ordinates of points a and b are... Then the matrix element mesh intersects with the fracture surface; the intersection point o The calculation formula is:

[0050]

[0051] 2) Among all the found intersection points, consider polygons formed by 3-8 intersection points, and then number the intersection points (i.e., the vertices of the polygon) in counter-clockwise order. (See appendix) Figure 4 This example uses a polygon formed by six intersection points to illustrate the rule for numbering these intersection points. The coordinates of the center point of the polygon formed by the six intersection points in the diagram are shown below. , , This indicates that the line connecting the intersection point to the center point forms an angle with the positive X-axis. The angle at intersection point 1 can be expressed as:

[0052] The same applies to the other points. All the calculated included angles are numbered from 1, 2, ..., n in ascending order, which is the counterclockwise order.

[0053] 3) Divide the generated intersection point numbers into a mesh of interstitial cells in sequence. (See attached image) Figure 5 The hexagonal mesh is divided into fracture element meshes in the order 1→2→3→4, 4→5→6→1, forming two fracture element meshes. Finally, all the divided fracture meshes are rotated back to the original coordinate system. This method yields a fracture finite element mesh with 4345 nodes and 5037 elements. See appendix for details. Figure 7 .

[0054] 4) Check the quality of the finite element mesh. The finite element mesh in this project example contains 4345 nodes and 5037 elements, most of which are quadrilateral meshes. The values ​​of the Jacobian determinants of the elements are all greater than 0 during isoparametric transformation, indicating that the mesh quality meets the requirements. Reasonably select the mesh size and element size to balance the accuracy and efficiency of finite element calculation, thereby improving work efficiency and ensuring the correctness and reliability of calculation data.

[0055] A geological model of an actual engineering project was created using traditional methods based on geological structural characteristics. The main adverse geological bodies considered in the 3D geological modeling process were large-scale faults (F14, F16, and F17), inter-layer fault zones (C5, C4, C3, C3-1, and C2) below the normal water level of 825m, and penetrating intra-layer fault zones LS331 and RS331 developed within columnar jointed rock strata P2β33. The model also included an arch dam structure, a left bank abutment, an enlarged foundation, and a seepage barrier. The 3D mesh was generated using the finite element analysis preprocessing software HyperMesh, and converted to a mesh data file (.nas) format acceptable to the COMSOL Multiphysics finite element software using its built-in interface. The model mesh primarily used tetrahedral elements, containing 231,890 nodes and 1,357,927 elements, as shown in the attached diagram. Figure 8 and Figure 9 As shown.

[0056] Compared with traditional methods, the method of this invention significantly reduces the number of meshes during mesh generation, achieving high-precision and high-efficiency generation of embedded fracture meshes, and establishing a mapping relationship between fracture meshes and three-dimensional matrix meshes.

[0057] The above embodiments are merely preferred technical solutions of the present invention and should not be considered as limitations on the present invention. The embodiments and features described in these embodiments can be arbitrarily combined without conflict. The scope of protection of the present invention should be limited to the technical solutions described in the claims, including equivalent substitutions of the technical features described in the claims. That is, equivalent substitutions and improvements within this scope are also within the scope of protection of the present invention.

Claims

1. A method for generating embedded finite element meshes for a three-dimensional fracture-matrix model, characterized in that: It includes the following steps: 1) Given a matrix 3D finite element mesh, determine the intersection between the fracture surface and the generated matrix 3D finite element mesh based on the spatial location and orientation of the fracture, and find all the intersection points between them; 2) Determine the coordinates of the centroid of the polygon formed by all intersection points, and number all intersection points in counterclockwise order; 3) Form a fracture element mesh according to the numbering order, and record the three-dimensional matrix mesh element numbers that intersect with the fracture elements; 4) Inspect the quality of the embedded crack finite element mesh; Step 1) includes: using the coordinate rotation method to find the intersection point between the fracture surface and the three-dimensional finite element mesh of the matrix, the specific steps of which are as follows: a. Let the plane containing any two intersecting edges of the fracture plane be the XOY plane, and let the endpoints A of the two edges be... B C Point A is the origin of the coordinate system, A→B is the positive X-axis, and the direction perpendicular to the XOY plane is the Z-axis. A new three-dimensional coordinate system is established. b. The coordinates of the finite element mesh nodes are changed from the original coordinate system. To the new coordinate system The transformation formula is: ； In the formula, the rotation matrix R can be expressed as: ； Matrix elements , , Let them be three pairwise orthogonal unit vectors, specifically represented as follows: set up It is the normal vector of the crack plane along the positive Z-axis, starting from point A. Recorded as , Recorded as Therefore , Then we have: ； because , can be , Represented as: ； ； Will Normalization yields: ； ； Finally, find the direction vector of the Y-axis: ； ； Then Normalization yields: ； If let middle, for , for , for ;but ； c. Determine whether any of the 12 edges of the matrix element mesh intersect the fracture surface. In the new coordinate system, iterate through the 12 edges of the matrix mesh element. If the first edge intersects the fracture surface... The coordinates of the two endpoints of the edge are respectively and ,and Then the line segment must intersect the crack plane. Let the coordinates of the intersection point be... ,in , The calculation formula is: ； ； Step 2) includes: if there is There are several intersection points, and the coordinates of these intersection points in the new coordinate system are as follows: , , ..., The coordinates of the center point of the enclosed planar figure are The calculation formula is: ； ； The distance from any intersection point to the center point is The x-coordinate of the intersection point relative to the center point is Then the angle between the vector from the center point to the intersection point and the positive x-axis is: ； If the ordinate of the intersection point is less than the ordinate of the center point Then the angle between it and the positive x-axis should be . Number the intersection points of all vectors with the positive x-axis in ascending order from 1, 2, ..., n, which is the counterclockwise order.

2. The embedded finite element mesh generation method for a three-dimensional fracture-matrix model according to claim 1, characterized in that: Step 3) includes: generating a fracture element mesh by sequentially generating the intersection point numbers. The maximum number of intersection points is 8. If the number of intersection points is less than 2, an element mesh cannot be formed. Considering the formation of a polygonal fracture element mesh by 3-8 intersection points, there are a total of 6 cases. The number of the three-dimensional matrix mesh element intersecting with the fracture element is recorded. If there are 3 intersection points, then the intersection points within the matrix grid cells will be combined into 1 fracture element grid, with the element numbering order being: 1→2→3→1; If the number of intersections is 4, then the intersections within the matrix grid cells will be combined into one fracture element grid, with the element numbering order being: 1→2→3→4; If the number of intersections is 5, then the intersections within the matrix mesh element will be combined into 2 fracture element meshes, with the element numbering order as follows: 1→2→3→4, 4→5→1→1; If the number of intersections is 6, then the intersections within the matrix grid cells will be combined into 2 fracture element grids, with the element numbering order as follows: 1→2→3→4, 4→5→6→1; If the number of intersections is 7, then the intersections within the matrix grid cells will be combined into 3 fracture element grids, with the element numbering order as follows: 1→2→3→4, 4→5→6→1, 6→7→1→1; If the number of intersections is 8, then the intersections within the matrix grid cells will be combined into 3 fracture element grids, with the element numbering order as follows: 1→2→3→4, 4→5→6→1, 6→7→8→1.

3. The embedded finite element mesh generation method for a three-dimensional fracture-matrix model according to claim 1 or 2, characterized in that: Step 3) further includes: calculating the coordinate values ​​of the fracture mesh nodes in the original coordinate system. The calculation formula is: 。 4. The embedded finite element mesh generation method for a three-dimensional fracture-matrix model according to claim 1, characterized in that: Step 4) includes: checking the quality of the finite element mesh. The crack element mesh generated by the program may have individual elements with a Jacobian determinant of less than 0. The mesh quality is improved by manually adjusting the mesh shape.

Images