Construction, evaluation and seepage flow simulation method of three-dimensional coarse discrete crack network

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 of fractured rock masses, and specifically relates to a method for constructing, evaluating, and simulating seepage flow of a three-dimensional rough discrete fracture network.

Background Art

[0002] Natural rock masses form numerous joints and fractures through long-term continuous actions and complex geological structure changes, greatly changing the mechanical properties and seepage flow characteristics of the rock masses and affecting the construction safety of underground engineering. Therefore, creating a real three-dimensional discrete-rough fracture network based on the occurrence of fractures in the field has become the focus of disaster prevention and mitigation research for underground engineering.

[0003] In recent years, with the development of computer technology and numerical calculation methods, the application of numerical calculation to the problems of rock mechanics and rock engineering has been increasing. From the perspective of the structure and form of rock masses, based on the statistics of a large amount of field measured data of rock fractures, the method of discrete fracture network modeling (DFN) using probability theory and mathematical statistics theory has been applied and verified in many actual projects. The methods and processes for generating 2D discrete fracture network models, visualization, and solving permeability coefficients have been widely studied. However, in the real world, fracture networks exist as 3D fractures, so there are errors in 2D discrete fracture network models. Therefore, Priest, Wang Enzhi, Zhang Guoqiang, Xie Jing, etc. have modeled 3D fracture networks and generated 3D fracture networks using disc models and square models by generating fracture parameters (density, length, orientation, dip angle, fracture aperture). In addition, numerical calculation research on the seepage flow characteristics of 3D fracture networks is developing. Scholars have studied seepage flow characteristics using mainstream numerical simulation software based on the CFD method, such as Ansys Fluent, OpenFOAM, Comsol Multiphysics, etc. Such commercial software is highly developed, has complete theories, and can basically meet the needs of most simulations. On the other hand, by using programming platforms such as C++, Java, Matlab, etc., available 3D seepage flow simulation software has been independently developed and studied. Currently, the Galerkin method, 3D Unified pipe-network Method (UPM), and VOF method are widely used.

[0004] However, the current research has the following three problems: (1) Three-dimensional cracks are often assumed to be smooth flat cracks, but actual cracks are rough, and the roughness of the cracks greatly changes the water permeability of a single network crack. Current research cannot consider the roughness of the crack surface. (2) Cracks are randomly distributed in a three-dimensional cracked rock mass, with large variability in distribution and a tendency to have a large number. When comparing two crack simulation methods, if the number of cracks is too large, business software often cannot directly model them, and there are many difficulties in assigning crack attribute values and setting boundary conditions, requiring desired secondary development. Independently developed simulation programs often have problems such as high calculation speed and high convergence, but low development level and narrow application range. (3) The research on the process of the seepage flow characteristics of a crack network includes the generation of the crack network, the calculation of connectivity, and the solution of the permeability coefficient. For the complex real crack network world, solving only a single crack network model at a time is inefficient, lacking a large-scale automatic solution method.

[0005] Therefore, currently, in order to facilitate the modeling, evaluation, and simulation of three-dimensional rough crack networks and contribute to actual water-related projects, a complete technical method for modeling three-dimensional discrete rough cracked rock masses based on the occurrence of rock strata, evaluating the connectivity of cracks, and simulating seepage flow is required.

Summary of the Invention

[0006] The main object of the present invention is to provide a method for constructing, evaluating, and simulating seepage flow of a three-dimensional rough discrete crack network in order to solve the problem lacking in the modeling method of a three-dimensional rough crack network, the evaluation of crack connectivity, and the method for simulating seepage flow in the prior art.

[0007] To achieve the above object, the present invention specifically provides a method for constructing, evaluating, and simulating seepage flow of a three-dimensional rough discrete crack network, including the following steps.

[0008] Step S1: Obtain a probability distribution model of the geometric parameters of the crack network based on field investigation data, boring data, geophysical exploration data, and other measured data.

[0009] Step S2: Use the programming software Matlab to randomly sample the parameters of the center Point_c of the crack disk, the crack dip angle Dip, the crack orientation Orientation, the radius Radius of the crack disk, and the crack aperture Aperture based on the Monte-Carlo method, and generate a data set of each parameter in the crack network.

[0010] Step S3: Use a random form of the Weierstrass function to generate height values at each coordinate of the upper and lower walls of the rough crack, and obtain a set of rough crack aperture widths for all cracks based on this.

[0011] Step S4: Read the data set of each parameter in the crack network obtained in Step S2, and evaluate the connectivity of the three-dimensional random crack network by the bounding box method.

[0012] Step S5: Based on the development platform Comsol with Matlab, perform coordinate system transformation of the crack disk and assignment of rough crack aperture values by the Euler angle and rotation matrix method, and create a numerical model of the three-dimensional discrete-rough crack network.

[0013] Step S6: Perform a simulation of the seepage flow of the three-dimensional crack network, and calculate the permeability coefficient of the three-dimensional crack 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 dip angle model, a crack orientation model, a crack disk radius model, and a crack aperture model.

[0015] Furthermore, each dataset in the crack network in step S2 is specifically as follows. Point_c={x 1 ,y 1 ,z 1 ;x 2 ,y 2 ,z 2 ;x 3 ,y 3 ,z 3 ;···;x n ,y n ,z n} Dip={α 1 ;α 2 ;α 3 ;···;α n} Orientation={β 1 ;β 2 ;β 3 ;···;β n} Radius={r 1 ;r 2 ;r 3 ;···;r n} Aperture={b 1 ;b 2 ;b 3 ;···;b n} n = L 長さ ·L 幅 ·L 高さ ·ρ where n is the number of cracks in the crack model, L 長さ , L 幅 and L 高さ represent the length, width, and height of the crack model respectively, x n , y n , z n represent the coordinates in the Cartesian coordinate system of the center point of the nth crack, α n is the inclination angle of the nth crack, β n is the orientation of the nth crack, r n is the radius of the nth crack disk, b n is the aperture 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 walls of the rough crack using a random-form Weierstrass function. The height of the upper crack wall is as follows.

Equation

[0018] Step S3.2: Read the crack aperture set Aperture = {b 1 ; b 2 ; b 3 ; ···; b n} in the crack network in step S2, and set Zupper(x i , y j )’ = Zupper(x i , y j ) + b 1 , that is, assign the initial aperture b i to the generated crack dataset Zupper(x j )’, and Zupper(x 1 , y i ) is the crack dataset without the initial aperture assigned. j ) is the crack dataset without the initial aperture assigned.

[0019] Step S3.3: Considering the crustal stress effect, input the stress and the normal stiffness to obtain the deformation amount of crack closure, and obtain the deformation amount of the actual rough crack opening width and the contact condition of the upper and lower wall surfaces.

Equation

[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 Equation (3).

Equation

[0021] Step S3.5: Repeat Steps S3.1 - S3.4 to obtain the set of rough crack - to - crack opening widths of all cracks.

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

[0023] Step S4.1: Initially judge the spatial relationship between the first crack and the second crack by the bounding box method, compare the straight - line distance L vb between the centers of the two boundary spheres and the radius of the spheres, analyze the overlap of the boundary spheres of the two cracks, and the geometric characteristics satisfy Equation (4).

Equation

[0024] Step S4.2: L vb > (a 1 + a 2 ) case, the normal vector n 1 = (l, m, n) is expressed as follows.

Equation

[0025] Step S4.3: The plane where the first cracked disk is located is expressed as follows. l 1 (x - x 1 ) + m 1 (y - y 1 ) + n 1 (z - z 1 ) = 0 (6) However, l 1 , m 1 , n 1 are the direction vectors of the first cracked disk.

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

Equation

[0027] Similarly, the plane and the boundary contour line where the second cracked disk is located can be expressed as follows respectively. l 2 (x - x 2 ) + m 2 (y - y 2 ) + n 2 (z - z 2 ) = 0 (8) (However, l 2 , m 2 , n 2(which is the direction vector of the second cracked disk)

Number

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

Number

[0029] Step S4.5: When the included angle θ = 0°, it indicates that the planes where the two cracked disks are located are parallel and there is no intersection line between the cracked disks. Otherwise, calculate whether the two cracked disks intersect with the intersection line respectively.

[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 cracked disk with the center of the disk as the origin of the plane global coordinate system and a radius of R 1 .

[0032] Step S5.2: Read the value b(x i , y j ) of the crack opening width between the newly generated upper and lower rough crack surfaces in Step S3, input the command into Comsol with Matlab, which is the development platform, and assign the value to the cracked disk so that the cracked disk has a rough crack opening width.

[0033] Step S5.3: Translate the origin of the global coordinate system parallel to the center point O of the disk, that is, the parallel translation directions along the x, y, and z directions of the origin of the global coordinate system are x 1 , y 1 , z 1 respectively.

[0034] Step S5.4: Introduce Euler's angle theorem, rotate the coordinate system by (180-α)° around the z-axis so that the x-axis of the global coordinate system is 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 - S5.4, iterate sequentially, and complete the modeling of the discrete - coarse crack network in three-dimensional space.

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

[0037] Step S6.1: Based on the soft platform Comsol with Matlab, call the governing equations and solver of the soft platform. The steady flow during the seepage flow simulation of the three-dimensional crack network is governed by Darcy's equation, specifically as follows.

[0038]

Equation

Equation

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

[0040] Step S6.3: Use the simulation software Comsol to mesh the model and solve the flow field and pressure field.

[0041] Step S6.4: Calculate the seepage flow velocity at the inlet of the crack network model, and obtain the steady-state volume flow rate Q under different crack aperture distributions by integrating the flow velocity at the inlet. Specifically, it is as follows. Q = Au (14) Where A = ∫dh, which is the cross-sectional area perpendicular to the fluid flow direction, and u is the velocity of the fluid passing through the cross-sectional area.

[0042] Step S6.5: Calculate the permeability coefficient K of the crack network model of the model according to Darcy's equation. Specifically, it is as follows.

Equation

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

[0044] Generate the parameters of the crack network based on the occurrence of rock cracks, and solve the problem of coordinate system conversion of the crack disk in Comsol, a business software platform, by using Euler angles and rotation matrices.

[0045] Realize the modeling and visualization of a three-dimensional "discrete - rough" crack network.

[0046] Automatically realize the process from the modeling of the crack network to the evaluation of the connectivity of the crack network and further to the seepage flow simulation in one go, reduce the cost of manual work, improve the calculation efficiency, and enable large-scale numerical simulation research.

[0047] To more clearly explain the embodiments for carrying out the present invention or the technical solutions in the prior art, the drawings necessary for the description of the embodiments for carrying out the invention or the prior art will be briefly described below. Obviously, the following drawings are some embodiments of the present invention, and those skilled in the art can obtain other drawings based on these drawings without creative efforts.

Brief Description of the Drawings

[0048]

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Embodiments for Carrying Out the Invention

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

[0050] The method for constructing, evaluating, and simulating seepage flow of the three-dimensional rough discrete fracture network shown in FIG. 1 specifically includes the following steps.

[0051] Step S1: Obtain a probability distribution model of the geometric parameters of the crack network as shown in Table 1 based on field investigation data, boring data, geophysical exploration data, and other measured data.

[0052]

Table 1

[0053] Step S2: Use the programming software Matlab to randomly sample the parameters of the center Point_c of the crack disk, the crack dip angle Dip, the crack orientation Orientation, the radius Radius of the crack disk, and the crack aperture Aperture based on the Monte-Carlo method, and generate a dataset of each parameter in the crack network as shown in Table 2.

[0054]

Table 2

[0055] Step S3: Use a random form of the Weierstrass function to generate height values at each coordinate of the upper and lower walls of the rough crack, and based on this, obtain a set of rough crack aperture widths for all cracks.

[0056] Step S4: Read the dataset of each parameter in the crack network obtained in Step S2, and evaluate the connectivity of the three-dimensional random crack network by the bounding box method.

[0057] Step S5: Based on the development platform Comsol with Matlab, perform coordinate system transformation of the crack disk and assignment of rough crack aperture values by the Euler angle and rotation matrix method, and create a numerical model of the three-dimensional discrete-rough crack network.

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

[0059] Specifically, conduct on-site investigation and research, and obtain the probability distribution model of the crack network parameters based on the measured data such as field investigation data, boring data, and geophysical exploration data, and determine the parameter values of each probability distribution model. In the present invention, the Baecher crack disk model is adopted, and the probability distribution models of the geometric parameters of the crack network 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 aperture width model.

[0060] The crack disk center model conforms to the Poisson distribution, and its formula is as follows.

Number

[0061] The crack inclination angle model and the crack orientation model conform to the Fisher distribution, and their formulas are as follows.

Number

[0062] The crack disk radius model conforms to the power distribution, and its formula is as follows. f(x)=x -α-1

[0063] The crack aperture width model conforms to the uniform distribution function, and its formula is as follows.

Number

[0064] Specifically, each data set in the crack network in step S2 is specifically as follows. Point_c={x 1 ,y 1 ,z1 ; x 2 , y 2 , z 2 ; x 3 , y 3 , z 3 ; ···; x n , y n , z n} Dip = {α 1 ; α 2 ; α 3 ; ···; α n} Orientation = {β 1 ; β 2 ; β 3 ; ···; β n} Radius = {r 1 ; r 2 ; r 3 ; ···; r n} Aperture = {b 1 ; b 2 ; b 3 ; ···; b n} = L 長さ · L 幅 · L 高さ · ρ In the formula, n is the number of cracks in the crack model, L 長さ , L 幅 and L 高さ represent the length, width, and height of the crack model respectively, x n , y n , z n represent the coordinates in the Cartesian coordinate system of the center point of the nth crack, α n is the inclination angle of the nth crack, β n is the orientation of the nth crack, r n is the radius of the nth crack disc, b n is the aperture width of the nth crack, and ρ is the crack density of the crack model.

[0065] Specifically, in the real world, cracks are rough, and the crack space (i.e., crack aperture) between the rough upper and lower walls of the crack affects the seepage flow characteristics within the crack. Currently, most research on rough cracks has been conducted on single cracks, and there is relatively little research on the modeling and seepage flow characteristics of rough cracks in a three-dimensional discrete crack network. To generate a three-dimensional discrete rough crack network, first, it is necessary to generate a dataset of the aperture between each rough crack for each crack. Taking the generation process of a single rough crack as an example, step S3 specifically includes the following steps.

[0066] Step S3.1: In order to make the generated crack aperture more in line with the actual situation, it is first necessary to generate the rough crack walls above and below a single crack. The roughness of the crack surface can be represented by the fractal dimension. Using the random form of the Weierstrass function, the height value at each coordinate of the upper and lower walls of the rough crack is generated, and the height of the upper crack wall is as follows.

[0067]

Equation

[0068] In the process of generating the upper and lower walls of the rough crack, while keeping other parameters constant, change the random number seed of C N to generate datasets of the upper and lower walls of different rough cracks Zupper(x i , y j ) and Zlower(x i , y j) can be generated.

[0069] Step S3.2: Read the crack opening width set Aperture = {b 1 ; b 2 ; b 3 ; ···; b n} in the crack network in step S2, and set Zupper(x i , y j )’ = Zupper(x i , y j ) + b 1 , that is, assign the initial opening width b i , y j )’ to the generated crack dataset Zupper(x 1 ), and Zupper(x i , y j ) is the crack dataset without the initial opening width assigned.

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

Number

[0071] Step S3.4: Lower (move down) Zupper(x i , y j ) by Δb f . That is, Zupper(x i , y j ) = Zupper(x i , y j ) - Δb f . Recalculate b(x i , y j ), that is, b(x i , y j ) = Zupper(xi , y j ) - Zlower(x i , y j ), and b(x i , y j ) = Zupper(x i , y j ) - Zupper(x i , y j ). If it is less than or equal to 0, it indicates that the upper and lower crack walls are in contact at this position, the crack closes at this position, and the fluid cannot pass through. Then, set b(x i , y j ) = 0. Update the data to obtain a single set of rough crack aperture data. Calculate the average crack aperture of the rough crack after the upper and lower cracks are displaced and partially closed by Equation (3).

Equation

[0072] Step S3.5: Repeat Steps S3.1 - S3.4 to obtain the set of rough crack apertures for all cracks.

[0073] Specifically, read the parameter set of each crack network and evaluate the connectivity of the three - dimensional random crack network. In the seepage flow process of the cracked rock mass, since isolated cracks cannot substantially affect the fluid movement process, the fluid flow in the crack network mainly depends on the connectivity of the rock mass crack network. The connectivity of the cracks is represented by the number of crack intersections and lines, and the connectivity index CI can be defined as follows.

Equation

[0074] Step S4 specifically includes the following steps.

[0075] Step S4.1: To calculate the intersection frequency of the cracked disks, the spatial relationship between the first crack and the second crack is initially determined by the bounding box method, and the linear distance L between the centers of the two boundary spheres vb and the radius of the sphere are compared, and the overlap of the boundary spheres of the two cracks is analyzed. As shown in FIGS. 2 and 3, the geometric features satisfy Equation (4).

Equation

[0076] Step S4.2: When L vb >(a 1 + a 2 ), it is necessary to further calculate and determine the positional relationship between the two cracks. According to the occurrence state of the cracked disk (which can be converted according to the inclination angle α and azimuth β, the radian value of the inclination angle, and the radian value of the azimuth), the normal vector n 1 =(l, m, n) of the cracked disk is expressed as follows.

Equation

[0077] Step S4.3: The plane where the first cracked disk is located is expressed as follows. l 1 (x - x 1 ) + m 1 (y - y 1 ) + n 1 (z - z 1 ) = 0 (6) However, l 1, m 1 , n 1 is the direction vector of the first cracked disk.

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

Number

[0079] Similarly, the plane where the second cracked disk is located and the boundary contour line can be expressed as follows, respectively. l 2 (x - x 2 ) + m 2 (y - y 2 ) + n 2 (z - z 2 ) = 0 (8) (where l 2 , m 2 , n 2 is the direction vector of the second cracked disk)

Number

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

Number

[0081] Step S4.5: When the included angle θ = 0, it indicates that the planes where the two cracked disks are located are parallel and there is no intersection line between the cracked disks. Otherwise, use equations (5) to (10) to calculate whether the two cracked disks intersect the intersection line respectively. If both of the two cracked disks intersect the intersection line, it is necessary to analyze and judge the positional relationship of the cracks based on the four spatial intersection points between the two cracked disks and the intersection line. As shown in FIGS. 4 and 5, the second crack being included in or intersecting the first crack means that the two cracked disks are in an intersecting relationship. At the same time, the spatial equation of the intersection line and the intersection length can be obtained using the four coordinates. FIG. 6 shows that the two cracks are separated, and the two cracked 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 cracked disk as an example, assume the parameters of the first cracked disk are O(x 1 , y 1 , z 1 ), inclination angle α 1 and azimuth β 1 , and disk radius R 1 respectively.

[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, set the origin of the plane global coordinate system as the center of the cracked disk, create a cracked disk with a radius of R 1 , and realize the visualization of the three-dimensional discrete crack network model as shown in FIG. 7.

[0085] Step S5.2: As shown in FIG. 8, the value b(x i , yj ) is read, commands are input into Comsol with Matlab which is a development platform, and its value is assigned to the cracked disk, enabling the assignment and visualization of the crack opening width value of the rough crack so that the cracked disk has a rough crack opening width.

[0086] Step S5.3: Translate the origin of the global coordinate system to the center point O of the disk, that is, the translation directions along the three directions of x, y, and z of the origin of the global coordinate system are x 1 , y 1 , z 1 respectively.

[0087] Step S5.4: Introduce Euler's angle theorem, rotate the coordinate system by (180-α)° around the z-axis so that the x-axis of the global coordinate system is 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 cracked disk model in the Cartesian coordinate system.

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

[0089] Specifically, step S6 specifically includes the following steps. Step S6.1: Based on Comsol with Matlab which is a software platform, call the governing equation and solver of the software platform. The steady flow during the seepage flow simulation of the three-dimensional crack network is governed by the Darcy equation, specifically as follows.

[0090]

Equation

Equation

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

[0092] Step S6.3: Use Comsol, which is simulation software, to mesh the model and solve for the flow field and pressure field.

[0093] Step S6.4: Calculate the seepage flow velocity at the inlet of the crack network model, and obtain the steady - state volume flow rate Q under different crack aperture distributions by integrating the flow velocity at the inlet. Specifically, it is as follows. Q = Au (14) In the formula, A = ∫dh, which is the cross - sectional area perpendicular to the fluid flow direction, and u is the velocity of the fluid passing through the cross - sectional area.

[0094] Step S6.5: Calculate the permeability coefficient K of the crack network model of the model according to Darcy's equation. Specifically, it is as follows.

Equation

[0095] The above is the modeling, connectivity evaluation and seepage flow solution process of a single three - dimensional rough crack network. By creating a script program and continuously reading the test plan, a parameter set of the crack network is generated and a program loop is performed.

[0096] Of course, the above description does not limit the present invention. The present invention is not limited to the above examples, and any 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 of obtaining a probability distribution model of geometric parameters of a fracture network based on 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 using programming software Matlab 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 by a random type Weierstrass function, and obtaining a set of rough crack opening widths for all cracks based on the height values; Step S4 of reading the data set of each parameter in the crack network obtained in step S2 and evaluating the connectivity of the three-dimensional random crack network by the bounding box method; Step S5: based on the development platform Comsol with Matlab, transform the coordinate system of the crack disk and assign the coarse crack opening width value by Euler angle and rotation matrix method to create a numerical model of a three-dimensional discrete-coarse crack network; A 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 in Matlab, and the origin of the plane global coordinate system is set as the center of the crack disk, and the radius is R １ Step S5.1 of creating a crack disk of The crack opening width value b(x ｉ , y ｊ ) and inputting commands into a development platform Comsol with Matlab to assign the value to the crack disk so that the crack disk has a rough crack opening width; The origin of the global coordinate system is translated to the center point O of the disk. In other words, 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 the Euler angle theorem, rotates the coordinate system around the z-axis by (180-α) degrees so that the x-axis of the global coordinate system is rotated to the orientation direction of the crack, and rotates the coordinate system around the y-axis by β so that the y-axis of the global coordinate system overlaps with the orientation line of the crack, thereby creating a crack disk model in the Cartesian coordinate system; and step S5.5 of repeating steps S5.1 to S5.4 and iteratively completing modeling of a discrete-coarse crack network in three-dimensional space; Specifically, step S6 is A step of invoking the governing equations and solvers of the software platform Comsol with Matlab, where the steady flow in the seepage flow simulation of the three-dimensional crack network is governed by the Darcy equation, specifically: [0027] (In the formula, d ｆ is the crack opening width, ρ is the crack density, μ is the dynamic viscosity, Q ｍ is the crack flow rate, u is the velocity, p is the pressure, k ｆ is the hydraulic conductivity, [0028] 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 at the inlet and outlet boundaries of the model, respectively; Step S6.3 of meshing the model and solving the flow and pressure fields using simulation software Comsol; 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) where A=∫dh, the cross-sectional area perpendicular to the direction of fluid flow, and u, the velocity of the fluid through the cross-sectional area; and Calculate the hydraulic conductivity K of the crack network model according to the Darcy equation, specifically: [0029] and step S6.5, where A is the cross-sectional area perpendicular to the fluid flow direction, L is the length of the crack 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 crack network as described in claim 1, characterized in that the probability distribution models of the geometric parameters of the crack network in step S1 include a crack disk circle center model, a crack inclination angle model, a crack orientation model, a crack disk radius model and a crack 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, 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 disk, and b ｎ The method for constructing, evaluating and simulating seepage flow in a three-dimensional coarse discrete crack network as described in claim 1, characterized in that: n is the opening width of the nth crack, and ρ is the crack density of the crack model.

4. Specifically, step S3 is The height values ​​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 [0030] (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 uniform distribution in [0, 2π], 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 ｊ )' with 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, inputting stress and normal stiffness to obtain the deformation amount of the crack closure, and obtaining the actual rough deformation amount of the crack opening width and the contact state of the upper and lower wall surfaces; [0031] (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 according to equation (3); [0032] (where M and N represent the number of nodes in the x and y directions, respectively) and step S3.5 of repeating steps S3.1 to S3.4 to obtain a set of coarse inter-crack opening widths for all cracks.

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 step S4.1, comparing the radius of the spheres, analyzing the overlap of the boundary spheres of the two cracks, and determining whether the geometric characteristics 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 crack disks, respectively.) L ｖｂ >(a １ +a ２ )in the case of, [0034] (l, m, n are directional vectors, β is orientation, and α is inclination angle) and the normal vector n of the crack disk １ = (l, m, n) in step S4.2; l １ (x-x １ )+m １ (y-y １ )+n １ (z-z １ )=0 (6) (However, １ , m １ , n １ is the direction vector of the first crack disk), which represents the plane on which the first crack disk is located, [Equation 35] Let denote 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 [0036] step S4.3, which can respectively represent the plane on which the second crack disk lies and the boundary contour line; Step S4.4: when the two boundary spheres overlap, calculate the included angle θ between the planes on which the two crack disks are located according to 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 of the crack discs; otherwise, calculating whether the two crack discs intersect with the intersection line or not.

Images