Methods for constructing, evaluating and simulating seepage flow in three-dimensional coarse-grained discrete fracture networks

Abstract

To solve the problem lacking in a modeling method of a three-dimensional coarse discrete crack network, crack connectivity evaluation, and a seepage flow simulation method in a conventional technique.SOLUTION: A construction, evaluation, and seepage flow simulation method of a three-dimensional coarse discrete crack network includes the following steps: acquiring a probability distribution model of geometric parameters of a crack network from outdoor investigation data, boring data, geophysical exploration data, and other actual measurement data, and generating data sets of each parameter in the crack network; generating height values of each coordinate of upper and lower wall surfaces of coarse cracks, and acquiring a set of coarse inter-crack opening widths of all the cracks based on the height values; reading out the data sets of each parameter and evaluating the connectivity of a three-dimensional random crack network using a bounding box method; generating a numerical model of the three-dimensional discrete-coarse crack network, performing a numerical simulation, and determining the permeability coefficient of the model.SELECTED DRAWING: Figure 1

Description

[Technical Field]

[0001] The present invention relates to the technical field of seepage flow calculation in fractured rock masses, and more particularly to a method for constructing, evaluating and simulating seepage flow in a three-dimensional coarse-grained discrete fracture network. [Background technology]

[0002] Natural rock masses undergo diagenesis and complex geological structural changes over a long period of time, resulting in the formation of numerous joints and fractures, which significantly change the mechanical properties and seepage characteristics of the rock mass and affect the safety of underground construction work. Therefore, creating an actual 3D discrete-coarse fracture network based on the occurrence of fractures in the field has become a key focus of research into disaster prevention and mitigation for underground construction work.

[0003] In recent years, advances in computer technology and numerical calculation methods have led to an increasing application of numerical calculations to problems in rock mechanics and rock engineering. Based on the statistical analysis of a large amount of field-measured data on rock fractures, discrete fracture network modeling (DFN) methods using probability theory and mathematical statistics have been applied and verified in many actual construction projects. Methods and processes for generating 2D discrete fracture network models, visualizing them, and solving for hydraulic conductivity have been widely studied. However, in the real world, fracture networks exist as 3D fractures, which means that 2D discrete fracture network models are prone to errors. For this reason, Priest, Wang Enzhi, Zhang Guoqiang, and Xie Jing et al. have modeled 3D fracture networks, generating fracture parameters (density, length, orientation, dip angle, and crack opening width) and using disk and square models to generate 3D fracture networks. Numerical research on the seepage flow characteristics of 3D fracture networks has also progressed. Scholars have used mainstream CFD-based numerical simulation software such as Ansys Fluent, OpenFOAM, and Comsol Multiphysics to study seepage flow characteristics. These commercial software programs are highly developed, have complete theories, and can basically meet most simulation needs. On the other hand, there has been independent development and research into available 3D seepage flow simulation software using programming platforms such as C++, Java, and Matlab. Currently, the Galerkin method, 3D Unified Pipe-Network Method (UPM), and VOF method are widely used.

[0004] However, current research has three limitations. (1) While 3D fractures are often assumed to be smooth, flat fractures, in reality, they are rough. The roughness of the fracture surfaces significantly affects the permeability of a single network of fractures, and current research fails to consider the roughness of the fracture surfaces. (2) Cracks are randomly distributed in 3D fractured rock masses, tending to be characterized by large variability in distribution and large numbers. Comparing the two fracture simulation methods, it is often difficult to directly model large numbers of fractures using commercial software. Assigning fracture attribute values and setting boundary conditions presents significant challenges, necessitating secondary development. Independently developed simulation programs offer fast calculation speeds and high convergence, but often suffer from low development levels and limited applicability. (3) Research into the process of seepage flow characteristics in fracture networks involves generating fracture networks, calculating connectivity, and solving for hydraulic conductivity. Given the complex nature of real fracture networks, solving a single fracture network model at a time is inefficient, and large-scale automated solutions are lacking.

[0005] Therefore, a complete technical method for modeling 3D discrete coarse-fractured rock masses based on rock formation occurrence, evaluating fracture connectivity, and simulating seepage flow is currently needed to facilitate the modeling, evaluation, and simulation of 3D coarse-fractured rock networks and contribute to actual water-related projects. Summary of the Invention

[0006] The main objective of the present invention is to provide a method for constructing, evaluating, and simulating seepage flow of a three-dimensional coarse discrete fracture network in order to solve the problem of the lack of methods for modeling three-dimensional coarse fracture networks, evaluating fracture connectivity, and simulating seepage flow in the prior art.

[0007] To achieve the above objectives, the present invention provides a method for constructing, evaluating and simulating seepage flow of a three-dimensional coarse discrete crack network, which specifically includes the following steps:

[0008] Step S1: Obtain a probability distribution model of the geometric parameters of the fracture network from field survey data, drilling data, geophysical survey data and other measured data.

[0009] Step S2: Using the programming software Matlab, the parameters of the crack disk center Point_c, crack inclination angle Dip, crack orientation Orientation, crack disk radius Radius, and crack opening width Aperture are randomly sampled based on the Monte-Carlo method to generate a data set for each parameter in the crack network.

[0010] Step S3: Generate height values at each coordinate of the upper and lower wall surfaces of the rough cracks using a random Weierstrass function, and obtain a set of rough crack opening widths for all cracks based on the height values.

[0011] Step S4: The data set of each parameter in the crack network obtained in step S2 is read, and the connectivity of the three-dimensional random crack network is evaluated using the bounding box method.

[0012] Step S5: Based on the development platform Comsol with Matlab, the coordinate system of the cracked disk is transformed and the coarse crack opening width value is assigned by the Euler angle and rotation matrix method, and a numerical model of the three-dimensional discrete-coarse crack network is constructed.

[0013] Step S6: Simulate the seepage flow of the three-dimensional fracture network and calculate the permeability coefficient of the three-dimensional fracture network model.

[0014] Furthermore, the probability distribution model of the geometric parameters of the crack network in step S1 includes a crack disk center model, a crack inclination angle model, a crack orientation model, a crack disk radius model, and a crack opening width model.

[0015] Furthermore, the specific data sets in the crack network in step S2 are as follows: Point_c={x1,y1,z1;x2,y2,z2;x3,y3,z3;···;x n ,y n ,z n} Dip={α1;α2;α3;···;α n} Orientation={β1;β2;β3;···;β n} Radius={r1;r2;r3;···;r n} Aperture={b1;b2;b3;···;b n} n=L 長さ L 幅 L 高さ ρ where n is the number of cracks in the crack model and L 長さ , L 幅 and L 高さ are the length, width, and height of the crack model, respectively, and x n ,y n ,z n represents the coordinate of the center point of the nth crack in the Cartesian coordinate system, and α n is the inclination angle of the nth crack, and β n is the orientation of the nth crack, and r n is the radius of the nth cracked disc, and b n is the opening width of the nth crack, and ρ is the crack density of the crack model.

[0016] Furthermore, step S3 specifically includes the following steps.

[0017] Step S3.1: Generate height values at each coordinate of the upper and lower wall surfaces of the rough crack by the random Weierstrass function, and the height of the upper crack wall surface is as follows:

number

[0018] Step S3.2: Set the crack opening widths in the crack network in step S2. Aperture = {b1;b2;b3; ;b n} and Zupper(x i ,y j )'=Zupper(x i ,y j )+b1, that is, the generated crack data set Zupper(x i ,y j )' is assigned the initial opening width b1, and Zupper(x i ,y j ) is a dataset of cracks for which no initial opening width is assigned.

[0019] Step S3.3: Taking into account the crustal stress action, the deformation amount of the crack closure is obtained by inputting the stress and normal stiffness, and the actual deformation amount of the rough crack opening width and the contact condition of the upper and lower wall surfaces are obtained.

number

[0020] Step S3.4: Calculate the average crack opening width of the rough crack after the upper and lower cracks are displaced and partially closed by Eq. (3).

number

[0021] Step S3.5: Repeat steps S3.1 to S3.4 to obtain a rough set of inter-crack opening widths for all cracks.

[0022] Furthermore, step S4 specifically includes the following steps.

[0023] Step S4.1: The spatial relationship between the first and second cracks is initially determined by the bounding box method, and the linear distance L between the centers of the two boundary spheres is calculated. vb and the radius of the sphere are compared, and the overlap of the boundary spheres of the two cracks is analyzed, and the geometric characteristics satisfy Eq. (4).

number

[0024] Step S4.2:L vb > (a1 + a2), the normal vector n1 = (l, m, n) of the cracked disk is expressed as follows:

number

[0025] Step S4.3: The plane on which the first cracked disc is located is expressed as follows: l1(x-x1)+m1(y-y1)+n1(z-z1)=0 (6) where l1, m1, and n1 are the direction vectors of the first crack disk.

[0026] The boundary contour of the first cracked disk is expressed as follows:

number

[0027] Similarly, the plane on which the second cracked disk lies and the boundary contour can be expressed as follows, respectively: l2(x-x2)+m2(y-y2)+n2(z-z2)=0 (8) (where l2, m2, and n2 are the direction vectors of the second cracked disk)

number

[0028] Step S4.4: When the two boundary spheres overlap, calculate the included angle θ between the planes on which the two cracked disks are located by Equation (10).

number

[0029] Step S4.5: If the included angle θ=0°, the planes on which the two crack discs are located are parallel, indicating that there is no intersection line between the crack discs; otherwise, calculate whether the two crack discs intersect with the intersection line.

[0030] Furthermore, step S5 specifically includes the following steps.

[0031] Step S5.1: Read and load the geometric parameter set of the crack network in step S2 in Matlab, and create a crack disk with radius R1, with the origin of the planar global coordinate system as the center of the crack disk.

[0032] Step S5.2: The crack opening width b(x) between the upper and lower rough crack faces newly generated in step S3 is calculated. i,y j ) and input commands into the development platform Comsol with Matlab to assign the value to the crack disc so that the crack disc has a rough crack opening width.

[0033] Step S5.3: The origin of the global coordinate system is translated to the center point O of the disk, that is, the translation directions along the three directions x, y, and z of the origin of the global coordinate system are defined as x1, y1, and z1, respectively.

[0034] Step S5.4: Introduce the Euler angle theorem, rotate the coordinate system by (180-α)° around the z-axis so that the x-axis of the global coordinate system is rotated to be in the crack orientation direction, and rotate the coordinate system by β around the y-axis so that the y-axis of the global coordinate system overlaps with the crack orientation line, and create a crack disk model in the Cartesian coordinate system.

[0035] Step S5.5: Repeat steps S5.1 to S5.4, and iterate to complete the modeling of the discrete-coarse crack network in the three-dimensional space.

[0036] Furthermore, step S6 specifically includes the following steps.

[0037] Step S6.1: Based on the software platform Comsol with Matlab, the governing equations and solver of the software platform are invoked. The steady flow in the seepage flow simulation of the 3D fracture network is governed by the Darcy equation, specifically as follows:

[0038]

number

number

[0039] Step S6.2: Set the boundary conditions, with the x direction being the fluid flow direction, and set as pressure boundary conditions at the inlet and outlet boundaries of the model, respectively.

[0040] Step S6.3: The model is meshed using the simulation software Comsol, and the flow field and pressure field are solved.

[0041] Step S6.4: Calculate the seepage flow velocity at the inlet of the fracture network model, and obtain the steady-state volumetric flow velocity Q under different fracture opening width distributions by integrating the flow velocity at the inlet, specifically as follows: Q = Au (14) where A = ∫dh, the cross-sectional area perpendicular to the direction of fluid flow, and u is the velocity of the fluid passing through the cross-sectional area.

[0042] Step S6.5: Calculate the hydraulic conductivity K of the fracture network model according to Darcy's equation, specifically as follows:

number

[0043] The present invention has the following beneficial effects.

[0044] The parameters of the fracture network are generated based on the occurrence of rock fractures, and the problem of the transformation of the coordinate system of the fracture disk in the business software platform Comsol is solved by using Euler angles and rotation matrices.

[0045] Achieve modeling and visualization of 3D "discrete-coarse" crack networks.

[0046] The process from fracture network modeling to fracture network connectivity evaluation and seepage flow simulation can be automated in one go, reducing the cost of manual work, improving computational efficiency, and enabling large-scale numerical simulation research.

[0047] In order to more clearly describe the embodiments for implementing the present invention or the technical solutions in the prior art, the following will briefly describe the drawings that need to be used in the description of the embodiments for implementing the present invention or the prior art. Obviously, the drawings below are only some embodiments of the present invention, and those skilled in the art can obtain other drawings based on these drawings without any creative efforts. [Brief explanation of the drawings]

[0048] [Figure 1] 1 shows a flowchart of a method for constructing, evaluating, and simulating seepage flow in a three-dimensional coarse discrete crack network according to the present invention. [Figure 2] A schematic diagram of two boundary spheres separated is shown. [Figure 3] A schematic diagram of two intersecting boundary spheres is shown. [Figure 4] 1 shows a schematic diagram of cracks contained within each other. [Figure 5] A schematic diagram of intersecting cracks is shown. [Figure 6] A schematic diagram of a crack separating is shown. [Figure 7] A diagram of the 3D discrete crack network model created by step S5.1 is shown. [Figure 8] Schematic diagram of the assignment of crack opening width values for rough cracks obtained by step S5.2 is shown. DETAILED DESCRIPTION OF THE INVENTION

[0049] The technical solutions of the present invention will be described clearly and completely below with reference to the drawings, and it is to be understood that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments, and other embodiments that can be obtained by those skilled in the art based on the embodiments of the present invention without any creative efforts shall all fall within the protection scope of the present invention.

[0050] The method for constructing, evaluating, and simulating the seepage flow of a three-dimensional coarse-grained discrete crack network shown in Figure 1 specifically includes the following steps.

[0051] Step S1: Obtain a probability distribution model of the geometric parameters of the fracture network as shown in Table 1 using field survey data, drilling data, geophysical survey data, and other measured data.

[0052] [Table 1]

[0053] Step S2: Using the programming software Matlab, the parameters of the crack disk center Point_c, crack inclination angle Dip, crack orientation Orientation, crack disk radius Radius, and crack opening width Aperture are randomly sampled based on the Monte-Carlo method, and a data set of each parameter in the crack network is generated as shown in Table 2.

[0054] [Table 2]

[0055] Step S3: Generate height values at each coordinate of the upper and lower wall surfaces of the rough cracks using a random Weierstrass function, and obtain a set of rough crack opening widths for all cracks based on the height values.

[0056] Step S4: The data set of each parameter in the crack network obtained in step S2 is read, and the connectivity of the three-dimensional random crack network is evaluated using the bounding box method.

[0057] Step S5: Based on the development platform Comsol with Matlab, the coordinate system of the cracked disk is transformed and the coarse crack opening width value is assigned by the Euler angle and rotation matrix method, and a numerical model of the three-dimensional discrete-coarse crack network is constructed.

[0058] Step S6: Simulate the seepage flow of the three-dimensional fracture network and calculate the permeability coefficient of the three-dimensional fracture network model.

[0059] Specifically, field surveys and research are conducted to obtain probability distribution models of the crack network parameters based on actual measurement data such as field survey data, drilling data, and geophysical survey data, and the parameter values of each probability distribution model are determined. In this invention, the Baecher crack disk model is adopted, and the probability distribution models of the crack network geometric parameters in step S1 include a crack disk center model, a crack inclination angle model, a crack orientation model, a crack disk radius model, and a crack opening width model.

[0060] The cracked disk center model fits the Poisson distribution, and its formula is as follows:

number

[0061] The crack inclination angle model and crack orientation model are fitted to the Fisher distribution, and the formula is as follows:

number

[0062] The crack disk radius model fits a power law distribution, and its formula is: f(x)=x -α-1

[0063] The crack opening width model is fitted with a uniform distribution function, the formula of which is:

number

[0064] Specifically, the data sets in the crack network in step S2 are as follows: Point_c={x1,y1,z1;x2,y2,z2;x3,y3,z3;···;x n ,y n ,z n} Dip={α1;α2;α3;···;α n} Orientation={β1;β2;β3;···;β n} Radius={r1;r2;r3;···;r n} Aperture={b1;b2;b3;···;b n} =L 長さ L 幅 L 高さ ρ where n is the number of cracks in the crack model and L 長さ , L 幅 and L 高さ are the length, width, and height of the crack model, respectively, and x n ,y n ,z n represents the coordinate of the center point of the nth crack in the Cartesian coordinate system, and α n is the inclination angle of the nth crack, and β n is the orientation of the nth crack, and r n is the radius of the nth cracked disc, and b n is the opening width of the nth crack, and ρ is the crack density of the crack model.

[0065] Specifically, in the real world, cracks are coarse, and the crack spacing between the upper and lower coarse walls (i.e., the crack opening width) affects the seepage characteristics within the crack. Currently, most research on coarse cracks has been conducted on single cracks, and there has been relatively little research on the modeling and seepage characteristics of coarse cracks within a 3D discrete crack network. To generate a 3D discrete coarse crack network, it is first necessary to generate a data set of inter-coarse crack opening widths for each crack. Taking the generation process of a single coarse crack as an example, step S3 specifically includes the following steps:

[0066] Step S3.1: To make the generated crack opening width more realistic, we first need to generate the upper and lower rough crack walls of a single crack. The roughness of the crack surface can be expressed as a fractal dimension, and the height values at each coordinate of the upper and lower walls of the rough crack are generated using a random Weierstrass function. The height of the upper crack wall is as follows:

[0067]

number

[0068] In the process of generating the rough crack upper and lower walls, C N The random number seed of Zupper(x i ,y j ) and Zlower(x i ,y j) can be generated.

[0069] Step S3.2: Set the crack opening widths in the crack network in step S2. Aperture = {b1;b2;b3; ;b n} and Zupper(x i ,y j )'=Zupper(x i ,y j )+b1, that is, the generated crack data set Zupper(x i ,y j )' is assigned the initial opening width b1, and Zupper(x i ,y j ) is a dataset of cracks for which no initial opening width is assigned.

[0070] Step S3.3: Taking into account the crustal stress action, the deformation amount of the crack closure is obtained by inputting the stress and normal stiffness, and the actual deformation amount of the rough crack opening width and the contact condition of the upper and lower wall surfaces are obtained.

number

[0071] Step S3.4: Zupper(x i ,y j ) to Δb f Lower (move down) by Zupper(x i ,y j )=Zupper(x i ,y j )-Δb f b(x i ,y j ) is calculated again, that is, b(x i ,y j )=Zupper(x i ,y j )-Zlower(x i ,y j), and b(x i ,y j )=Zupper(x i ,y j )-Zupper(x i ,y j ) ≦ 0, the upper and lower crack walls come into contact, the crack closes at this position, and the fluid cannot pass through. i ,y j ) = 0. The data is updated to obtain a single coarse crack opening width data set. The average crack opening width of the coarse crack after the upper and lower cracks have been displaced and partially closed is calculated using equation (3).

number

[0072] Step S3.5: Repeat steps S3.1 to S3.4 to obtain a rough set of inter-crack opening widths for all cracks.

[0073] Specifically, the parameter set for each fracture network is read and the connectivity of the 3D random fracture network is evaluated. In the seepage flow process of fractured rock, isolated fractures cannot substantially affect the fluid movement process, so the fluid flow in the fracture network mainly depends on the connectivity of the rock fracture network. The connectivity of fractures is expressed by the number of fracture intersections and lines, and the connectivity index CI can be defined as follows:

number

[0074] Step S4 specifically includes the following steps:

[0075] Step S4.1: To calculate the frequency of cracked disk intersections, the spatial relationship between the first and second cracks is initially determined by the bounding box method, and the linear distance L between the centers of the two bounding spheres is calculated. vb By comparing the radius of the sphere and the boundary sphere of the two cracks, the overlap of the boundary spheres of the two cracks is analyzed, and the geometric characteristics satisfy Equation (4), as shown in Figures 2 and 3.

number

[0076] Step S4.2:L vb If >(a1+a2), the relative positions of the two cracks need to be further calculated and determined. According to the occurrence of the crack disk (which can be converted according to the dip angle α and the orientation β, the dip angle in radians, and the orientation in radians), the normal vector n1=(l,m,n) of the crack disk can be expressed as follows:

number

[0077] Step S4.3: The plane on which the first cracked disc is located is expressed as follows: l1(x-x1)+m1(y-y1)+n1(z-z1)=0 (6) where l1, m1, and n1 are the direction vectors of the first crack disk.

[0078] The boundary contour of the first cracked disk is expressed as follows:

number

[0079] Similarly, the plane on which the second cracked disk lies and the boundary contour can be expressed as follows, respectively: l2(x-x2)+m2(y-y2)+n2(z-z2)=0 (8) (where l2, m2, and n2 are the direction vectors of the second cracked disk)

number

[0080] Step S4.4: When the two boundary spheres overlap, calculate the included angle θ between the planes on which the two cracked disks are located by Equation (10).

number

[0081] Step S4.5: If the included angle θ = 0, the planes on which the two crack disks are located are parallel, indicating that there is no intersection line between the crack disks. Otherwise, calculate whether the two crack disks intersect with the intersection line using equations (5) to (10). If both crack disks intersect with the intersection line, the positional relationship of the cracks must be analyzed and determined based on the four spatial intersection points between the two crack disks and the intersection line. As shown in Figures 4 and 5, if the second crack is contained in the first crack or intersects with the first crack, it means that the two crack disks are in an intersection relationship. At the same time, the spatial equation of the intersection line and the intersection length can be calculated using the four coordinates. Figure 6 shows that the two cracks are separated and the two crack disks have the same intersection line but do not intersect. Here, A 11 and A 12 are the first and second intersection points between crack A and the intersection line, respectively, and B 11 and B 12 are the first and second intersection points between crack B and the intersection line, respectively.

[0082] Specifically, taking the generation process of the first crack disk as an example, the parameters of the first crack disk are assumed to be O(x1, y1, z1), inclination angle α1, orientation β1, and disk radius R1.

[0083] Step S5 specifically includes the following steps:

[0084] Step S5.1: Read and load the geometric parameter set of the crack network in step S2 in Matlab. The origin of the planar global coordinate system is set as the center of the crack disk, and a crack disk with a radius of R1 is created. As shown in Figure 7, the visualization of the 3D discrete crack network model is realized.

[0085] Step S5.2: As shown in Figure 8, the crack opening width value b(x) between the upper and lower rough crack faces newly generated in step S3 is calculated. i ,y j ) and input commands into the development platform Comsol with Matlab to assign the value to the crack disc, so that the crack disc has a coarse crack opening width, thereby realizing the assignment and visualization of the crack opening width value of the coarse crack.

[0086] Step S5.3: The origin of the global coordinate system is translated to the center point O of the disk, that is, the translation directions along the three directions x, y, and z of the origin of the global coordinate system are defined as x1, y1, and z1, respectively.

[0087] Step S5.4: Introduce the Euler angle theorem, rotate the coordinate system by (180-α)° around the z-axis so that the x-axis of the global coordinate system is rotated to be in the crack orientation direction, and rotate the coordinate system by β around the y-axis so that the y-axis of the global coordinate system overlaps with the crack orientation line, and create a crack disk model in the Cartesian coordinate system.

[0088] Step S5.5: Repeat steps S5.1 to S5.4, and iterate to complete the modeling of the discrete-coarse crack network in the three-dimensional space.

[0089] Specifically, step S6 specifically includes the following steps: Step S6.1: Based on the software platform Comsol with Matlab, the governing equations and solver of the software platform are invoked. The steady flow in the seepage flow simulation of the 3D fracture network is governed by the Darcy equation, specifically as follows:

[0090]

number

number

[0091] Step S6.2: Set the boundary conditions, with the x direction being the fluid flow direction, and set as pressure boundary conditions at the inlet and outlet boundaries of the model, respectively.

[0092] Step S6.3: The model is meshed using the simulation software Comsol, and the flow field and pressure field are solved.

[0093] Step S6.4: Calculate the seepage flow velocity at the inlet of the fracture network model, and obtain the steady-state volumetric flow velocity Q under different fracture opening width distributions by integrating the flow velocity at the inlet, specifically as follows: Q = Au (14) where A = ∫dh, the cross-sectional area perpendicular to the direction of fluid flow, and u is the velocity of the fluid passing through the cross-sectional area.

[0094] Step S6.5: Calculate the hydraulic conductivity K of the fracture network model according to Darcy's equation, specifically as follows:

number

[0095] The above is the process of modeling a single 3D coarse fracture network, evaluating its connectivity, and solving for seepage flow. A script program is created to generate a set of fracture network parameters by continuously reading the test plan, and the program loops.

[0096] Of course, the above description does not limit the present invention, and the present invention is not limited to the above examples, and any variations, modifications, additions or substitutions made by those skilled in the art within the substantial scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. specifically, Step S1: Obtaining a probability distribution model of geometric parameters of a fracture network from field survey data, drilling data, geophysical survey data, and other actual measurement data; Step S2: randomly sampling the parameters of the crack disk center Point_c, the crack inclination angle Dip, the crack orientation Orientation, the crack disk radius Radius, and the crack opening width Aperture based on the Monte-Carlo method to generate a data set of each parameter in the crack network; Step S3: generating height values at each coordinate of the upper and lower wall surfaces of the rough cracks using a random Weierstrass function, and obtaining a set of rough crack opening widths for all cracks based on the height values; Step S4: reading the data set of each parameter in the crack network acquired in step S2 and evaluating the connectivity of the three-dimensional random crack network using the bounding box method; Step S5: transforming the coordinate system of the crack disk and assigning coarse crack opening width values by Euler angles and rotation matrix method to create a numerical model of the 3D discrete-coarse crack network; and step S6 of simulating the seepage flow of the three-dimensional fracture network and calculating the hydraulic conductivity of the three-dimensional fracture network model; Specifically, step S5 is The geometric parameter set of the crack network in step S2 is read and loaded, and the origin of the plane global coordinate system is the center of the crack disk, and the radius is R １ Step S5.1 of creating a crack disc of The crack opening width b(x) between the upper and lower rough crack surfaces newly generated in step S3 ｉ , y ｊ ) and assigning the value to the crack disc so that the crack disc has a coarse crack opening width; The origin of the global coordinate system is translated to the center point O of the disk. That is, the translation directions of the origin of the global coordinate system along the three directions of x, y, and z are respectively defined as x １ , y １ , z １ Step S5.3: Step S5.4 introduces Euler angle theorem, rotates the coordinate system by (180-α)° around the z-axis so that the x-axis of the global coordinate system is rotated to the crack orientation direction, and rotates the coordinate system by β around the y-axis so that the y-axis of the global coordinate system overlaps with the crack orientation line, thereby creating a crack disk model in the Cartesian coordinate system; Step S5.5 of repeating steps S5.1 to S5.4 and iteratively completing the modeling of the discrete-coarse crack network in three-dimensional space; Specifically, step S6 is Invoking the governing equations and solver of the soft platform, the steady flow in the seepage flow simulation of the three-dimensional fracture network is governed by the Darcy equation, specifically: [0000] (In the formula, d ｆ is the crack opening width, ρ is the crack density, μ is the kinematic viscosity, Q ｍ is the crack flow rate, u is the velocity, p is the pressure, k ｆ is the permeability, [0000] is the gradient operator of calculus, which represents differentiation in different directions); Step S6.2 of setting boundary conditions, with the x-direction being the fluid flow direction and pressure boundary conditions being set at the inlet and outlet boundaries of the model, respectively; step S6.3 of meshing the model and solving the flow and pressure fields; The seepage flow velocity at the inlet of the crack network model is calculated, and the steady-state volumetric flow velocity Q under different crack opening width distributions is obtained by integrating the flow velocity at the inlet. Specifically, Q = Au (14) step S6.4, where A=∫dh, the cross-sectional area perpendicular to the direction of fluid flow, and u, the velocity of the fluid passing through the cross-sectional area; Calculate the hydraulic conductivity K of the crack network model according to Darcy's equation, specifically: [0000] and step S6.5, where A is the cross-sectional area perpendicular to the fluid flow direction, L is the length of the fracture network in the flow direction, and Δp is the pressure difference between the inlet and outlet boundaries of the model.

2. The method for constructing, evaluating and simulating seepage flow of a three-dimensional coarse discrete fracture network, as described in claim 1, characterized in that the probability distribution model of the geometric parameters of the fracture network in step S1 includes a fracture disk center model, a fracture inclination angle model, a fracture orientation model, a fracture disk radius model, and a fracture opening width model.

3. Specifically, each data set in the crack network in step S2 is Point_c={x １ ,y １ ,z １ ;x ２ ,y ２ ,z ２ ;x ３ ,y ３ ,z ３ ;・・・;x ｎ ,y ｎ ,z ｎ } Dip = {a} １ ;a ２ ;a ３ ;···;a ｎ } Orｰentation={b １ ;b ２ ;b ３ ;・・・;b ｎ } Radius={r １ ;r ２ ;r ３ ;・・・;r ｎ } Aperture={b １ ;b ２ ;b ３ ;・・・;b ｎ } n=L 長さ ・L 幅 ・L 高さ ・ρ (where n is the number of cracks in the crack model, and L 長さ , L 幅 and L 高さ are the length, width, and height of the crack model, respectively, and x ｎ , y ｎ , z ｎ represents the coordinate of the center point of the nth crack in the Cartesian coordinate system, and α ｎ is the inclination angle of the nth crack, and β ｎ is the orientation of the nth crack, and r ｎ is the radius of the nth cracked disc, and b ｎ The method for constructing, evaluating and simulating seepage flow in a three-dimensional coarse discrete fracture network as described in claim 1, characterized in that: ρ is the opening width of the nth fracture, and ρ is the fracture density of the fracture model.

4. Specifically, step S3 is The height values at each coordinate of the upper and lower walls of the rough crack are generated by the random Weierstrass function, and the height of the upper crack wall is [Equation 30] (In the formula, Zupper ｉ、ｊ (x ｉ , y ｊ ) is the roughness of the crack wall (x ｉ , y ｊ ) the height of the plane coordinates at Ｎ are independent random numbers with standard normal distribution, N is the number of random numbers, D and λ are fractal variables, A Ｎ and B Ｎ are independent random numbers with a uniform distribution in [0, 2π], and i, j are the spatial frequency resolution of the generated crack network), and C Ｎ The random number seed of the crack is changed to obtain the data set Zupper (x ｉ , y ｊ ) and Zlower(x ｉ , y ｊ ) and In step S2, the crack opening width set in the crack network is Aperture = {b １ ;b ２ ;b ３ ;...;b ｎ } and Zupper(x ｉ , y ｊ )'=Zupper(x ｉ , y ｊ ) + b １ That is, the data set of cracks to be generated is Zupper (x ｉ , y ｊ )' to the initial opening width b １ and Zupper(x ｉ , y ｊ ) is a data set of cracks for which no initial opening width is assigned; and Step S3.3: Taking into account the crustal stress action, the deformation amount of the crack closure is obtained by inputting the stress and normal stiffness, and the deformation amount of the actual rough crack opening width and the contact state of the upper and lower wall surfaces are obtained. [Equation 31] (In the formula, Δb ｆ is the deformation of the crack closure, V ｍ is the deformation amount of the initial crack closure, σ ｎ is the normal stress, K ｎｉ is the initial normal stiffness) Step S3.4: calculating the average crack opening width of the coarse crack after the upper and lower cracks are displaced and partially closed by equation (3); [Equation 32] (where M and N represent the number of nodes in the x and y directions, respectively) The method for constructing, evaluating, and simulating seepage flow of a three-dimensional coarse discrete fracture network as described in claim 1, further comprising: step S3.5 of repeating steps S3.1 to S3.4 to obtain a set of coarse inter-crack opening widths for all fractures.

5. Specifically, step S4 is The spatial relationship between the first and second cracks is initially determined by the bounding box method, and the linear distance L between the centers of the two boundary spheres is calculated. ｖｂ and the radius of the sphere, and the boundary sphere of the two cracks Step S4.1: Analyzing the overlap and determining whether the geometric features satisfy equation (4); [Equation 33] (In the formula, L ｖｂ is the distance between the boundary spheres of the first and second cracks, O １ (x １ , y １ , z １ ), O ２ (x ２ , y ２ , z ２ ) are the global coordinates of the centers of the first and second cracked disks, respectively, and a １ , a ２ represent the radii of the first and second cracked disks, respectively). L ｖｂ >(a １ +a ２ )in the case of, [Equation 34] (l, m, n are direction vectors, β is orientation, α is inclination angle) and the normal vector n of the crack disk １ step S4.2 representing = (l, m, n); l １ (x-x １ )+m １ (y-y １ )+n １ (z-z １ )=0 (6) (However, １ , m １ , n １ is the direction vector of the first cracked disc), which represents the plane on which the first cracked disc is located, [Equation 35] represents the boundary contour of the first cracked disk, Similarly, l ２ (x-x ２ )+m ２ (y-y ２ )+n ２ (z-z ２ )=0 (8) (However, ２ , m ２ , n ２ is the direction vector of the second cracked disk), and [Equation 36] Step S4.3, where the plane on which the second cracked disc lies and the boundary contour line can be represented by Step S4.4: when the two boundary spheres overlap, calculate the included angle θ of the plane where the two crack disks are located by equation (10); [Equation 37] The method for constructing, evaluating and simulating a three-dimensional coarse discrete crack network as described in claim 1, characterized in that it includes step S4.5: if the included angle θ = 0, the planes on which the two crack discs are located are parallel, indicating that there is no intersection line between the crack discs; otherwise, it includes step S4.5: respectively calculating whether the two crack discs intersect with the intersection line.

Images