Method for geometric processing of three-dimensional embedded discrete fracture network

Abstract

This invention discloses a three-dimensional embedded discrete crack network geometric processing method applicable to cracks with arbitrary tilt angles, including: reading crack data for each crack; calculating the coordinates of the four vertices of the crack; calculating crack segments with tilt angles not equal to 90° and equal to 90° respectively; this invention solves the problem of difficult geometric processing of cracks with arbitrary tilt angles in three-dimensional EDFM models.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology for oil and gas field development, and in particular to a geometric processing method for three-dimensional embedded discrete fracture networks. Background Technology

[0002] Conventional numerical simulation methods, such as the discrete element method (DEM), finite element method (FEM), and finite difference method (FDDM), suffer from complex mesh generation and enormous computational costs when simulating complex crack networks. In three-dimensional models, complex crack networks may even be impossible to mesh, or the resulting meshes may be of poor quality and fail to meet the basic mesh requirements of numerical simulation. Therefore, researchers have proposed the Embedded Discrete Crack Network Model (EDFM). This model extends the conventional numerical simulation method's view of cracks and matrix as a single mesh, separating them into two independent meshes, thus partially solving the problem of preprocessing complex crack networks.

[0003] The most crucial aspect of EDFM models is identifying the segments formed after a complete crack is divided by the matrix mesh. Currently, most EDFM model preprocessing methods proposed in research focus on 2D EDFM models, with only a small portion addressing 3D EDFM. Furthermore, most 3D EDFM research deals with the special case of a 90° inclination angle (the angle between the crack and the X-axis is the azimuth, and the angle with the XY plane is the inclination angle), with few researchers proposing universal algorithms for arbitrary inclination angles. In the 2D case, the matrix is ​​a plane, and the crack is a line segment. The crack is divided into several line segments by the matrix mesh, and the segmented crack segments can be easily identified by finding the coordinates of the intersection points between the short crack lines and the mesh. However, in the 3D EDFM model, the crack is elevated from a line segment to a plane. When the crack is perpendicular to the matrix surface, regardless of the azimuth angle, the resulting shape after the crack is divided by the matrix mesh is rectangular or parallelogram-shaped, making it easy to identify. However, when the inclination angle is not 90°, depending on the azimuth angle, the crack plane is divided into polygons containing 3-6 vertices by the matrix mesh, significantly increasing the difficulty of identification. Summary of the Invention

[0004] To address the challenge of geometric processing of cracks at arbitrary tilt angles in 3D EDFM models, this invention proposes a 3D embedded discrete crack network geometric processing method, which is a general geometric processing approach.

[0005] To achieve the above objectives, the technical solution adopted by this invention is: a three-dimensional embedded discrete crack network geometric processing method, applicable to cracks tilted at any angle, specifically including the following steps:

[0006] Step 1: Read the crack data for each crack, which includes the crack center coordinates Pc(x,y,z), azimuth angle α, tilt angle β, crack length L, crack height H, and crack width W.

[0007] Step 2: Using the crack data, calculate the coordinates of the four vertices p0, p1, p2 and p3 of the crack after two rotations; the first rotation is a counterclockwise rotation of the azimuth angle α around the Z-axis, and the second rotation is a clockwise rotation of the azimuth angle α based on the coordinate system obtained in the first rotation.

[0008] Step 3: When the inclination angle β of the crack is not equal to 90°, project the crack in the XY direction to obtain the projected image PF, and calculate {Xmin,Xmax,Ymin,Ymax} in the crack coordinates.

[0009] Step 4: Establish a set of mesh numbers MC corresponding to the current crack. Iterate through the X and Y coordinates of the center of all meshes. If the coordinates are within the range of {Xmin, Xmax, Ymin, Ymax}, then put the current mesh number into the set MC.

[0010] Step 5: Divide the set MC into two subsets, inside and cross. Iterate through the coordinates of the four vertices of the grid contained in the set MC. inSideNumber is the number of vertices of the current grid in the PF graph. If inSideNumber is equal to 4, put the grid number into inside; if inSideNumber is greater than 0 and less than 4, put the grid number into cross; if inSideNumber is equal to 0, remove the grid number from MC.

[0011] Step 6: Establish a set FS of the split crack segments. For the mesh inside, directly combine the 4 vertices in clockwise order to form a crack segment and put it into the set FS.

[0012] Step 7: Divide the edge part into two parts: vertices and intermediate intersections. Use the convex hull algorithm to process the mesh at the edge and put the found convex polygons into the set FS to form crack segments.

[0013] Step 8: Calculate the plane equation of the plane where the crack is located using any three points from the vertex coordinates obtained in Step 2. Traverse the set FS and substitute the X and Y coordinates of each vertex coordinate of all crack segments into the plane equation to obtain the correct Z coordinate.

[0014] Step 9: When the inclination angle β of the crack is equal to 90°, the intersection coordinates of p0-p1 with the grid are directly calculated by interpolation to obtain all the X and Y coordinates after the crack is divided. The intersection coordinates of p0-p3 with the grid are calculated to obtain all the Z coordinates after the crack is divided. The X, Y and Z coordinates are combined in order to obtain all the crack segments after the crack is divided.

[0015] Step 10: Calculate the coordinates of the intersection point of each crack segment with the Z-axis. Calculate whether each edge of the crack segment intersects with the coordinates of the intersection point containing the Z-axis. If it does, then divide the crack segment into two cracks using the convex hull algorithm.

[0016] As a further improvement of the present invention, step 2 specifically includes the following steps:

[0017] Step 2.1: Calculate the coordinates of the points pe(x,y,z) and pb(x,y,z) on the crack plane that intersect the crack edge with the straight line line passing through the crack center coordinate Pc and parallel to the XY plane using the crack length L. Then calculate the coordinates pb_p(x,y,z) and pe_p(x,y,z) after rotating counterclockwise by α.

[0018] Step 2.2: Calculate the projected coordinates p0_p, p1_p, p2_p, and p3_p of the crack apex in the rotated coordinate system using the coordinates of pb_p(x,y,z) and pe_p(x,y,z) combined with the crack length L, crack height H, and tilt angle β.

[0019] Step S2.3: Rotate the coordinate system clockwise by α to restore it, and obtain the coordinates (x, y, z) of each vertex of the crack in the original coordinate system.

[0020] As a further improvement to the present invention, the formula for calculating the coordinates before and after rotation is as follows:

[0021] First counterclockwise rotation: p_p = p*Rz(θ); Second clockwise rotation: p = p_p*Rz(-θ);

[0022] The rotation coordinate matrix is:

[0023]

[0024] As a further improvement of the present invention, in step 7, the convex hull algorithm adopts the Graham Scan algorithm to find the convex hull containing the most points in the initial point set, and then sorts the points to construct the convex hull using the Graham Scan algorithm, and finally outputs the coordinate set of the sorted points.

[0025] As a further improvement of the present invention, step 7 specifically includes the following steps:

[0026] Step 7.1: For the crack vertex, first extract the coordinates of the intersection point P_f that is closest to P in the clockwise direction and the coordinates of the intersection point P_b that is closest to P in the counterclockwise direction, to form the initial point set PC;

[0027] Step 7.2: Traverse the cross set constructed in step 5, calculate the distance between the center of each grid in cross and P, and the grid with the closest distance is the grid where P is located, and extract the coordinates of the four vertices of the grid.

[0028] Step 7.3: By determining the coordinates of the four vertices of the grid and the coordinates of P_f and P_b, the number of vertices of the convex polygon can be preliminarily determined to be (3, 4, 5).

[0029] Step 7.4: For the intersection points between crack vertices, the two adjacent intersection points are P_f and P_b, and the crack center Pc is the corresponding P. The initial point set PC contains only P_f and P_b.

[0030] Step 7.5: Sort the points in PC using the Graham Scan algorithm, and finally output the vertex coordinates contained in the crack fragment, and add the crack fragment to the FS set;

[0031] Step 7.6: Repeat steps 7.1-7.5 to process the crack segment at the crack edge.

[0032] As a further improvement of the present invention, step 7.3 further includes: further determining which vertices in the mesh belong to the crack segment, as follows:

[0033] 3 vertices: The crack segment does not contain mesh vertices, so no judgment is needed. P, P_f, and P_b directly constitute the crack segment.

[0034] Four vertices: A crack segment contains one mesh vertex. Construct two judgment points Pj1(P_f.x,P_b.y) and Pj2(P_b.x,P_f.y). The vertex in the mesh with the same coordinates as Pj1 or Pj2 is the vertex of the crack segment. Add this vertex to PC.

[0035] Five vertices: The crack segment contains two mesh vertices. Construct two vectors v1(P-P_f) and v2(P_P_b), and calculate s1 = v1 x v2; traverse the two vertices P of the mesh clockwise. ci and P ci+1 Construct two vectors V3(PP) ci ) and V4 (PP ci+1 ), calculate S2 = v3 x v4, if s1*s2>0, then P ci and P ci+For the vertices of the crack fragment, add 2 vertices to PC;

[0036] As a further improvement of the present invention, the method for determining the mesh vertices in step 7.4 is as follows:

[0037] Three vertices: The crack segment contains one mesh vertex. Construct two judgment points Pj1(P_f.x,P_b.y) and Pj2(P_b.x,P_f.y). The vertex in the mesh with the same coordinates as Pj1 or Pj2 is the vertex of the crack segment. Add this vertex to PC.

[0038] Four vertices: The crack segment contains two mesh vertices. Construct two vectors v1(P-P_f) and v2(P_P_b), and calculate s1 = v1 x v2; traverse the two vertices P of the mesh clockwise. ci and P ci+1 Construct two vectors V3(PP) ci ) and V4 (PP ci+1 ), calculate S2 = v3 x v4, if s1*s2>0, then P ci and P ci+ For the vertices of the crack fragment, add 2 vertices to PC.

[0039] As a further improvement to the present invention, the following steps are also included:

[0040] Step 11: If there are multiple cracks, repeat steps 1-10 to complete the geometric preprocessing of all cracks contained in the complex crack network.

[0041] The beneficial effects of this invention are:

[0042] This invention determines the meshed pattern of a crack by calculating its projection onto a plane, eliminating the need to calculate the projections of crack segments onto three separate planes. It considers the possibility of 3, 4, and 5 convex hull vertices when the crack edge intersects the mesh and proposes corresponding judgment algorithms. The proposed method for constructing crack point sets and their subsets processes all cracks by traversing the X and Y coordinates of the base mesh only once, rather than requiring a complete traversal of the X, Y, and Z coordinates for each crack, resulting in higher computational efficiency. Attached Figure Description

[0043] Figure 1 This is a schematic diagram of the calculation and projection of the four vertices of the crack after reading the crack parameters in an embodiment of the present invention;

[0044] Figure 2 This is a schematic diagram of the point selection for {Xmin, Xmax, Ymin, Ymax} in an embodiment of the present invention;

[0045] Figure 3 This is a schematic diagram illustrating the secondary partitioning of the MC set into inside and cross sets in an embodiment of the present invention;

[0046] Figure 4 This is a schematic diagram illustrating three cases in which the vertices of the crack edge form a convex polygon in an embodiment of the present invention;

[0047] Figure 5 This is a schematic diagram illustrating two cases in which the intersection points between the vertices of the crack edge form a convex polygon in an embodiment of the present invention;

[0048] Figure 6 This is a schematic diagram of the crack segment indicated by the crack apex and the middle intersection point in an embodiment of the present invention;

[0049] Figure 7 This is a schematic diagram of the calculation results after the crack segment is projected onto the XY plane and cut by the mesh in the XY direction in an embodiment of the present invention;

[0050] Figure 8 This is a schematic diagram of the calculation results after the crack is restored by projection onto the XY plane and then cut by the XYZ direction in an embodiment of the present invention;

[0051] Figure 9 This is a schematic diagram of the final calculation results of the matrix and crack mesh in an embodiment of the present invention. Detailed Implementation

[0052] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0053] Example

[0054] A geometric processing method for a three-dimensional embedded discrete crack network of cracks with arbitrary angles includes the following steps:

[0055] The matrix mesh portion is generated based on the basic geometric information of the mesh:

[0056] The reservoir has a grid size of 50m*25m*5m, with a grid division of (20, 20, 20) in the XYZ directions, totaling 8000 grids. A fracture exists in the reservoir with the following parameters: Pc(10, 12.5, 2.5), azimuth 15°, dip 60°, fracture length L 10m, fracture height H 4m, and fracture width 0.1m (these parameters do not affect the calculation results).

[0057] S10: Read the crack center coordinates Pc(x,y,z), azimuth angle α, tilt angle β, crack length L, crack height H, and crack width W for each crack (this data structure is the common crack representation method for most crack simulation software).

[0058] S20: Based on the read data, the coordinates of the four vertices of the crack are first calculated (the crack is generally rectangular or elliptical; an ellipse can be approximated by dividing it into several rectangles). Because the crack has an azimuth angle, the coordinates of the crack vertices p0, p1, p2, and p3 are obtained by rotating the coordinate system twice. The first rotation is a counterclockwise rotation of the azimuth angle around the Z-axis, and the second rotation is a clockwise rotation of the azimuth angle based on the coordinate system obtained in the first rotation. Figure 1 (The direction of rotation is shown in the diagram). The specific calculation process is as follows:

[0059] Step S21: First, calculate the coordinates of the points pe(x,y,z) and pb(x,y,z) on the crack plane that intersect the crack edge with the straight line line passing through PC and parallel to the XY plane through L. Then calculate the coordinates pb_p(x,y,z) and pe_p(x,y,z) after rotating counterclockwise by α.

[0060] Step S22: At this point, the straight line is in the XZ plane and parallel to the X-axis. The crack edges p0-p3 and p1-p2 are parallel to the YZ plane. By combining the coordinates of pb_p(x,y,z) and pe_p(x,y,z) with L, H, and β, we can easily calculate the projected coordinates p0_p, p1_p, p2_p, and p3_p of the crack apex in the rotated coordinate system. After obtaining the projected coordinates ( Figure 1 The black hollow circle in the middle is the projection point;

[0061] Step S23: Rotate the coordinate system clockwise by α to restore it, obtaining the coordinates (x, y, z) of each vertex of the crack in the original coordinate system. The formulas for calculating the coordinates before and after rotation are:

[0062] (1) First counterclockwise rotation: p_p=p*Rz(θ);

[0063] (2) Second clockwise rotation; p = p - p * Rz(-θ)

[0064] The rotation coordinate matrix is:

[0065]

[0066] The calculated vertex coordinates are shown in Table 1.

[0067] Table 1

[0068] Vertex Number p0 p1 p2 p3 coordinate (5.429,10.239,0.767) (15.088,12.8281,0.767) (14.570,14.76,4.232) (4.911,12.171,4.232) Vertex Number p0_p p1_p p2_p p3_p coordinate (7.894,8.485,0.767) (17.894,8.485,0.767) (17.894,10.485,4.232) (7.894,10.485,4.232)

[0069] S30: First, handle the general case where β is not equal to 90°. Project the crack into the XY direction to obtain the projected image PF, and calculate {Xmin, Xmax, Ymin, Ymax} in the crack coordinates, such as... Figure 2As shown, the calculated actual results are {5.429, 7.894, 8.485, 4.232}.

[0070] S40: Establish a set of grid numbers MC corresponding to the current crack. Iterate through the X and Y coordinates of the grid centers of all grids. If the coordinates are within the range of {Xmin, Xmax, Ymin, Ymax}, then add the current grid number to the set MC. In this way, when multiple cracks coexist, all cracks can be processed by iterating through the grid center coordinates only once, instead of re-traversing the entire grid for each crack, thus improving algorithm efficiency. The MC set is shown as the shaded area of ​​the grid in Figure 2, and the calculated MC set is shown in Table 2.

[0071] Table 2

[0072]

[0073] S50: Further process the MC set for each crack, dividing the MC into two subsets: inside and cross. Iterate through the four vertex coordinates of the meshes contained in the MC. `inSideNumber` represents the number of vertices of the current mesh within the PF graph. If `inSideNumber` equals 4, add the mesh number to the inside subset; if `inSideNumber` is greater than 0 and less than 4, add the mesh number to the cross subset; if `inSideNumber` equals 0, remove the mesh number from the MC. The set after this secondary partitioning is as follows: Figure 3 As shown, the quadratic partitioning result is an empty set inside, and the cross set (the calculation result is the same as MC) is shown in Table 3:

[0074] Table 3

[0075]

[0076] S60: Establish a set FS of the split crack segments. For the mesh inside, its 4 vertices are directly one of the graphics after the PF is split by the crack. Directly assemble the 4 vertices into a crack segment in a clockwise order (the internal split part must be a rectangle or parallelogram) and put it into FS.

[0077] S70: After processing the inside part, the remaining part of the PF after being divided by the mesh is the edge of the crack. The edge part is further divided into two parts: vertices and intermediate intersections. By analyzing the shape of the crack edge after being divided by the mesh through an "exhaustive" approach, it is found that the shape divided at the edge must be a convex polygon (the PF must also be a convex polygon). The mesh at the edge is processed using the convex hull algorithm, and the crack segments composed of the found convex polygons are put into the FS.

[0078] The convex hull algorithm uses the Graham Scan algorithm. This algorithm takes a set of points as input, finds the convex hull containing the most points in the set, and finally outputs the coordinates of the sorted points. Because the number of vertices of the convex polygon formed after the crack edge is divided by a mesh is already determined when the mesh intersections are determined, the calculation of the convex polygon requires directly determining the point set and then using the Graham Scan algorithm to sort the points and construct the convex hull. The specific process for determining the initial point set is as follows:

[0079] Step S71: For the crack vertex, first extract the coordinates of the intersection point P_f that is closest to P in the clockwise direction and the coordinates of the intersection point P_b that is closest to P in the counterclockwise direction, to form the initial point set PC;

[0080] Step S72: Traverse the `cross` set constructed in step S40, calculate the distance between the center of each grid in `cross` and P, and the grid with the closest distance is the grid containing P. Extract the coordinates of the four vertices of this grid;

[0081] Step S73: By determining the coordinates of the four vertices of the mesh and the coordinates of P_f and P_b, the number of vertices of the convex polygon can be initially determined (3, 4, 5). Based on this, it is necessary to further determine which vertices in the mesh belong to the crack segment. The method for determining the number of vertices differs depending on the specific number of vertices, as follows (e.g., ...). Figure 4 As shown):

[0082] (1) 3 vertices: The crack segment does not contain mesh vertices, so no judgment is needed. P, P_f, and P_b directly constitute the crack segment.

[0083] (2) 4 vertices: The crack segment contains 1 mesh vertex. Construct two judgment points Pj1(P_f.x,P_b.y) and Pj2(P_b.x,P_f.y). The vertex in the mesh with the same coordinates as Pj1 or Pj2 is the vertex of the crack segment. Add the vertex to PC.

[0084] (3) 5 vertices: The crack segment contains 2 mesh vertices. Construct two vectors v1(P-P_f) and v2(P_P_b), and calculate s1 = v1 x v2. Traverse the two vertices P of the mesh clockwise. ci and P ci+1 Construct two vectors V3(PP) ci ) and V4 (PP ci+1 ), calculate S2 = v3 x v4, if s1*s2>0, then P ci and P ci+ For the vertices of the crack fragment, add 2 vertices to PC;

[0085] Step S74: For the intersection points between crack vertices, the two adjacent intersection points are P_f and P_b, and the crack center Pc is the corresponding P. The initial point set PC only contains P_f and P_b. The methods for determining the mesh vertices are as follows (e.g.) Figure 5 As shown):

[0086] (1) 3 vertices: The crack segment contains 1 mesh vertex. Construct two judgment points Pj1(P_f.x,P_b.y) and Pj2(P_b.x,P_f.y). The vertex in the mesh with the same coordinates as Pj1 or Pj2 is the vertex of the crack segment. Add the vertex to PC.

[0087] (2) Four vertices: The crack segment contains two mesh vertices. Construct two vectors v1(P-P_f) and v2(P_P_b), and calculate s1 = v1 x v2. Traverse the two vertices P of the mesh clockwise. ci and P ci+1 Construct two vectors V3(PP) ci ) and V4 (PP ci+1 ), calculate S2 = v3 x v4, if s1*s2>0, then P ci and P ci+ For the vertices of the crack fragment, add 2 vertices to PC;

[0088] Step S75: Sort the points in PC using the Graham Scan algorithm, and finally output the vertex coordinates contained in the crack fragment, and add the crack fragment to the FS set.

[0089] S80: Repeat steps S70-S75 to process the crack segment at the crack edge.

[0090] S90: Currently, all the crack segments found in FS are images projected onto the XY plane (only the X and Y coordinates are correct), and they still need to be restored. Calculate the plane equation of the plane where the crack is located using any three points from the vertex coordinates obtained in step S20. Traverse FS and substitute the X and Y coordinates from the vertex coordinates of all crack segments into the plane equation to obtain the correct Z coordinate.

[0091] S100: When β is 90°, the intersection coordinates of p0-p1 with the grid are directly calculated by interpolation to obtain all the X and Y coordinates after the crack is divided. The intersection coordinates of p0-p3 with the grid are calculated to obtain all the Z coordinates after the crack is divided. By combining the X, Y and Z coordinates in order, all the crack segments after division can be obtained.

[0092] S110: After processing the XY directions, calculate the coordinates of the intersection point of each crack segment with the Z axis, and calculate whether each edge of the crack segment intersects with the coordinates of the Z axis. If it does, then divide the crack segment into two cracks using the convex hull algorithm.

[0093] S120: If multiple cracks exist, repeat steps S10-S110 to complete the geometric preprocessing of all cracks contained in the complex crack network. The calculation is now complete, and the final matrix and crack mesh generation results are as follows: Figure 9 As shown.

[0094] The embodiments described above are merely illustrative of specific implementations of the present invention, and while the descriptions are detailed, they should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention.

Claims

1. A geometric processing method for a three-dimensional embedded discrete crack network, characterized in that, Applicable to cracks tilted at any angle, specifically including the following steps: Step 1: Read the crack data for each crack, which includes the crack center coordinates Pc(x,y,z), azimuth angle α, tilt angle β, crack length L, crack height H, and crack width W. Step 2: Using the crack data, calculate the coordinates of the four vertices p0, p1, p2 and p3 of the crack after two rotations; the first rotation is a counterclockwise rotation of the azimuth angle α around the Z-axis, and the second rotation is a clockwise rotation of the azimuth angle α based on the coordinate system obtained in the first rotation. Step 3: When the inclination angle β of the crack is not equal to 90°, project the crack in the XY direction to obtain the projected image PF, and calculate {Xmin,Xmax,Ymin,Ymax} in the crack coordinates. Step 4: Establish a set of mesh numbers MC corresponding to the current crack. Iterate through the X and Y coordinates of the center of all meshes. If the coordinates are within the range of {Xmin, Xmax, Ymin, Ymax}, then put the current mesh number into the set MC. Step 5: Divide the set MC into two subsets, inside and cross. Iterate through the coordinates of the four vertices of the grid contained in the set MC. inSideNumber is the number of vertices of the current grid in the PF graph. If inSideNumber is equal to 4, put the grid number into inside; if inSideNumber is greater than 0 and less than 4, put the grid number into cross; if inSideNumber is equal to 0, remove the grid number from MC. Step 6: Establish a set FS of the split crack segments. For the mesh inside, directly combine the 4 vertices in clockwise order to form a crack segment and put it into the set FS. Step 7: Divide the edge part into two parts: vertices and intermediate intersections. Use the convex hull algorithm to process the mesh at the edge and put the found convex polygons into the set FS to form crack segments. Step 8: Calculate the plane equation of the plane where the crack is located using any three points from the vertex coordinates obtained in Step 2. Traverse the set FS and substitute the X and Y coordinates of each vertex coordinate of all crack segments into the plane equation to obtain the correct Z coordinate. Step 9: When the inclination angle β of the crack is equal to 90°, the intersection coordinates of p0-p1 with the grid are directly calculated by interpolation to obtain all the X and Y coordinates after the crack is divided. The intersection coordinates of p0-p3 with the grid are calculated to obtain all the Z coordinates after the crack is divided. The X, Y and Z coordinates are combined in order to obtain all the crack segments after the crack is divided. Step 10: Calculate the coordinates of the intersection point of each crack segment with the Z-axis. Calculate whether each edge of the crack segment intersects with the coordinates of the intersection point containing the Z-axis. If it does, then divide the crack segment into two cracks using the convex hull algorithm.

2. The three-dimensional embedded discrete crack network geometric processing method according to claim 1, characterized in that, Step 2 specifically includes the following steps: Step 2.1: Calculate the coordinates of the points pe(x,y,z) and pb(x,y,z) on the crack plane that intersect the crack edge with the straight line line passing through the crack center coordinate Pc and parallel to the XY plane using the crack length L. Then calculate the coordinates pb_p(x,y,z) and pe_p(x,y,z) after rotating counterclockwise by α. Step 2.2: Calculate the projected coordinates p0_p, p1_p, p2_p, and p3_p of the crack apex in the rotated coordinate system using the coordinates of pb_p(x,y,z) and pe_p(x,y,z) combined with the crack length L, crack height H, and tilt angle β. Step S2.3: Rotate the coordinate system clockwise by α to restore it, and obtain the coordinates (x, y, z) of each vertex of the crack in the original coordinate system.

3. The three-dimensional embedded discrete crack network geometric processing method according to claim 2, characterized in that, The formulas for calculating the coordinates before and after rotation are: First counterclockwise rotation: p_p = p*Rz(θ); Second clockwise rotation: p = p_p*Rz(-θ); The rotation coordinate matrix is:

4. The three-dimensional embedded discrete crack network geometric processing method according to claim 1, 2, or 3, characterized in that, In step 7, the convex hull algorithm uses the Graham Scan algorithm to find the convex hull containing the most points in the initial point set. Then, the points are sorted and the convex hull is constructed using the Graham Scan algorithm. Finally, the coordinate set of the sorted points is output.

5. The three-dimensional embedded discrete crack network geometric processing method according to claim 4, characterized in that, Step 7 specifically includes the following steps: Step 7.1: For the crack vertex, first extract the coordinates of the intersection point P_f that is closest to P in the clockwise direction and the coordinates of the intersection point P_b that is closest to P in the counterclockwise direction, to form the initial point set PC; Step 7.2: Traverse the cross set constructed in step 5, calculate the distance between the center of each grid in cross and P, and the grid with the closest distance is the grid where P is located, and extract the coordinates of the four vertices of the grid. Step 7.3: By determining the coordinates of the four vertices of the grid and the coordinates of P_f and P_b, the number of vertices of the convex polygon can be preliminarily determined to be (3, 4, 5). Step 7.4: For the intersection points between crack vertices, the two adjacent intersection points are P_f and P_b, and the crack center Pc is the corresponding P. The initial point set PC contains only P_f and P_b. Step 7.5: Sort the points in PC using the Graham Scan algorithm, and finally output the vertex coordinates contained in the crack fragment, and add the crack fragment to the FS set; Step 7.6: Repeat steps 7.1-7.5 to process the crack segment at the crack edge.

6. The three-dimensional embedded discrete crack network geometric processing method according to claim 5, characterized in that, Step 7.3 also includes: further determining which vertices in the mesh belong to the crack fragment, as follows: 3 vertices: The crack segment does not contain mesh vertices, so no judgment is needed. P, P_f, and P_b directly constitute the crack segment. Four vertices: A crack segment contains one mesh vertex. Construct two judgment points Pj1(P_f.x,P_b.y) and Pj2(P_b.x,P_f.y). The vertex in the mesh with the same coordinates as Pj1 or Pj2 is the vertex of the crack segment. Add this vertex to PC. Five vertices: The crack segment contains two mesh vertices. Construct two vectors v1(P-P_f) and v2(P_P_b), and calculate s1 = v1 x v2; traverse the two vertices P of the mesh clockwise. ci and P ci+1 Construct two vectors V3(PP) ci ) and V4 (PP ci+1 ), calculate S2 = v3 x v4, if s1*s2>0, then P ci and P ci+ For the vertices of the crack fragment, add 2 vertices to PC.

7. The three-dimensional embedded discrete crack network geometric processing method according to claim 6, characterized in that, In step 7.4, the method for determining the mesh vertices is as follows: Three vertices: The crack segment contains one mesh vertex. Construct two judgment points Pj1(P_f.x,P_b.y) and Pj2(P_b.x,P_f.y). The vertex in the mesh with the same coordinates as Pj1 or Pj2 is the vertex of the crack segment. Add this vertex to PC. Four vertices: The crack segment contains two mesh vertices. Construct two vectors v1(P-P_f) and v2(P_P_b), and calculate s1 = v1 x v2; traverse the two vertices P of the mesh clockwise. ci and P ci+1 Construct two vectors V3(PP) ci ) and V4 (PP ci+1 ), calculate S2 = v3 x v4, if s1*s2>0, then P ci and P ci+ For the vertices of the crack fragment, add 2 vertices to PC.

8. The three-dimensional embedded discrete crack network geometric processing method according to claim 7, characterized in that, It also includes the following steps: Step 11: If there are multiple cracks, repeat steps 1-10 to complete the geometric preprocessing of all cracks contained in the complex crack network.

Images