A numerical simulation method for estimating equivalent seepage of three-dimensional fracture network based on topological properties

Abstract

The application discloses a numerical simulation method for estimating equivalent seepage of a three-dimensional fracture network based on topological properties, which is based on relevant data of a fracture network of a research area, first constructs a small three-dimensional fracture network sample set, then maps the sample fracture network into an equivalent two-dimensional undirected graph by using graph theory, solves equivalent permeability, and constitutes a training sample set; then the maximum information index method is used to mine main control factors of the equivalent permeability of the fracture network, and the control relationship between the fracture network parameters and the underground flow field is clarified; finally, a substitution model between the main control factors and the equivalent permeability is constructed based on a BPNN machine learning model, and conversion of a target discrete fracture network to an equivalent medium is realized. The application provides a simulation scheme with high prediction accuracy and low calculation cost reflecting topological connectivity for numerical simulation of an underground flow field of a fractured stratum at a site scale.

Description

Technical Field

[0001] This invention relates to a numerical simulation method for flow fields in fracture networks in underground strata, specifically a numerical simulation method for estimating equivalent seepage in a three-dimensional fracture network based on topological properties. Background Technology

[0002] Fracture networks are widely present in natural rock matrices, providing favorable paths for fluid flow and acting as water barriers, thus influencing the distribution of underground flow fields. Therefore, for well-consolidated bedrock, obtaining the equivalent permeability of fracture networks in underground strata, relative to their relatively low pore permeability, is crucial for seepage simulation. However, most popular fracture equivalent permeability models rely on observed geometric features while neglecting the topological structure of the fracture network. In fact, fracture distribution exhibits significant anisotropy, and the topological structure has a strong controlling effect on seepage. Discrete fracture models can accurately describe the connectivity within fracture networks and fluid flow within fractures, but their computational cost is high when the fracture network is complex. Equivalent permeable medium models are more computationally efficient, but they do not consider topological connectivity, thus significantly reducing model complexity, resulting in large deviations in calculation results. Summary of the Invention

[0003] Purpose of the invention: The purpose of this invention is to provide a numerical simulation method for estimating the equivalent seepage of a three-dimensional fracture network based on topological properties in fractured strata at the site scale, which has high prediction accuracy and low computational cost.

[0004] Technical solution: The numerical simulation method for estimating equivalent seepage flow in a three-dimensional fracture network based on topological properties, as described in this invention, includes:

[0005] (1) Collect basic data on fracture networks and basic data on hydrogeology, and construct a target discrete three-dimensional fracture network based on the basic data on fracture networks;

[0006] (2) Construct a conceptual model of the underground flow field of the fracture network based on the target discrete three-dimensional fracture network and basic hydrogeological data;

[0007] (3) Based on geological conditions, hydrogeological conditions and target area scale, determine the coverage range of fracture network parameters, and then use the Monte Carlo method to establish several small three-dimensional fracture networks by randomly generating fracture orientation to form an initial training sample set.

[0008] (4) Use graph theory to map the sample fracture network into an equivalent two-dimensional undirected graph, and solve for the equivalent permeability of the sample fracture network;

[0009] (5) The maximum information index method is used to mine the main control parameters of equivalent permeability in several fracture network parameters, construct a sample set with the main control parameters as input and equivalent permeability as output, and normalize the input and output variables.

[0010] (6) Based on the BPNN machine learning model, a prediction model of the equivalent permeability of a three-dimensional fracture network is constructed. The target discrete fracture network is gridded, and the three-dimensional matrix composed of the master control parameters is used as the input. The three-dimensional matrix composed of the equivalent permeability, i.e. the equivalent medium, is used as the output. The hyperparameters of the model are adjusted to achieve high-precision simulation of input and output. Finally, a numerical model for simulating the underground flow field of the fracture network is obtained.

[0011] Furthermore, in step (1), the basic data of the fracture network includes fracture orientation distribution characteristics, fracture diameter distribution characteristics, and fracture density distribution characteristics.

[0012] Furthermore, basic hydrogeological data includes structural geological maps and hydrogeological maps.

[0013] Furthermore, in step (3), the orientation parameters of the small three-dimensional fracture network are as follows:

[0014] The crack diameter follows a power-law distribution, expressed as:

[0015] f(l)=βl -α

[0016] Where f(l) is the probability density function of the crack diameter distribution; β is the value based on the minimum crack diameter l. min and maximum value l max The obtained normalization exponent; α is the power law exponent;

[0017] The crack diameter is expressed as:

[0018]

[0019] Where R is a random number;

[0020] The fracture orientation follows the Fisher distribution:

[0021]

[0022] Where θ is the dip angle; θ0 is the average dip angle; φ is the dip direction; φ0 is the average dip direction; κ is the dispersion parameter; and sinh is the hyperbolic sine function.

[0023] The crack aperture λ is between 0.1 mm and 1 mm, conforming to a uniform random distribution; it is assumed that the cracks in each crack network have the same aperture.

[0024] Further, in step (4), graph theory is used to map the sample crack network into an equivalent two-dimensional undirected graph, including:

[0025] Treating the cracks as nodes V and the intersections between cracks as edges E, the 3D crack network is transformed into a 2D undirected graph G(V, E) as follows:

[0026] The three-dimensional fracture network F consists of N fractures f, represented as:

[0027] F={f i}, i = 1, ..., N

[0028] Define a bijective mapping Φ:

[0029] Φ:f i →v i

[0030] If there are two cracks f i with f j intersect, Then there is an edge in E connecting the corresponding vertex:

[0031]

[0032] Among them, (v i v j )∈E represents vertex v i and v j The edge e between ij , side e ij weight w ij Set as vertex v i and v j Length of the line of intersection between them:

[0033] w ij =Lengthf i ∩f j

[0034] Similarly, considering the flow direction, we can regard the inlet plane x0 as the source node s, and the outlet plane x L Considering the target node t, when the crack intersects with the inlet and outlet boundaries:

[0035]

[0036]

[0037] Among them, e si For an edge between a source node and a non-source / sink node, e it The edges between non-source / sink nodes and the target node are defined; thus, the 3D fracture network is transformed into a 2D undirected graph.

[0038] When the 3D fracture network is sparse, the 2D undirected graph will be divided into several parts. Only the fracture sub-network corresponding to the connected subgraph containing source nodes and target nodes can effectively connect the inlet plane and the outlet plane and induce seepage. Therefore, the connected subgraph without source nodes and target nodes can be ignored. That is, the 2D undirected graph G can be pruned into the union of all connected subgraphs containing source nodes and target nodes except for edge nodes, i.e., a subset G' of the 2D undirected graph G. G' has a unique corresponding sub-network F' of the 3D fracture network F.

[0039] Furthermore, the equivalent permeability of the sample fracture network is solved, including: solving the equivalent permeability of the fracture network in the flow direction on the equivalent subgraph G' of the subnetwork F', based on the following assumptions:

[0040] i. The rock matrix is ​​incompressible and impermeable;

[0041] ii. The fluid is incompressible;

[0042] iii. The two planes of the crack are parallel and smooth;

[0043] iv. Fluid flow follows Darcy's law and exhibits laminar flow behavior.

[0044] Furthermore, assuming a two-dimensional undirected graph G has M nodes, where each node represents the intersection of fractures, and node i is adjacent to node j, according to the cube law, the total flow Q between node i and node j is... j Represented as:

[0045]

[0046] Where k is the equivalent permeability; g is the gravitational acceleration; μ is the hydrodynamic viscosity; I is the hydraulic gradient; λ ij b is the aperture between crack i and crack j; b is the crack width; H i For the upstream head, H j For downstream head, L ij Let be the Euclidean distance between node i and node j;

[0047] According to the law of conservation of mass, the flow rate at node i is expressed as:

[0048]

[0049] The regional head distribution of the sample fracture network is represented as follows:

[0050]

[0051] Where, ΔH ij The difference in water head between upstream and downstream; w ijLet w be the equivalent water conductivity between node i and node j, and also the weight of the edge connecting node i and node j. When there is no connection between node i and node j, w is... ij The value is 0;

[0052] From this, we can obtain the adjacency matrix A and the degree matrix D of the vertices of a two-dimensional undirected graph G, defined as follows:

[0053] The equivalent permeability k of the entire sample fracture network region is expressed as:

[0054]

[0055] Among them, L zone Q is the region length of the sample crack network. si It is the inflow, H s It's the inlet head, H t It's the water outlet head.

[0056] Furthermore, in step (5), the fracture network parameters include topological and geometric parameters, the topological parameters including:

[0057] Degree centrality represents the number of fractures that intersect with the fracture.

[0058] Eigenvector centrality is the iterative exponent of degree centrality, representing the importance of the fracture in the entire fracture network;

[0059] Clustering coefficient represents the diversity of pathways in a fracture network;

[0060] Betweenness centrality, which represents the number of the fracture in the shortest path of the network, plays an important role in controlling seepage.

[0061] Maximum flow is a feasible flow with the maximum capacity in a network with capacity, and it is an important topological feature of seepage flow.

[0062] Anisotropy index f A The expression is:

[0063]

[0064] Where N is the number of fractures, λ is the fracture aperture, i is the unit vector of the flow direction, and n is the fracture normal vector;

[0065] Geometric parameters include the number of fractures, fracture aperture, and fracture area P per unit volume. 32 and unit cell size.

[0066] Furthermore, in step (5), the maximum information index method is used to mine the master control parameters, including:

[0067] Suppose the dataset has Γ data points. Plot a scatter plot of the two attributes in a two-dimensional space, dividing it into a grid of size m×n along the horizontal and vertical axes, where m×n < B and B = Γ. 0.6 The grid with coordinates (x0, y0) contains γ data points, with a frequency P(x, y) = γ / Γ. Based on this, the maximum grid resolution of the normalized mutual information is calculated as the metric of MIC.

[0068]

[0069] Where X and Y are the two variables whose correlation needs to be calculated;

[0070] The MIC value ranges from |0 to 1|. When two variables are independent, the MIC approaches 0; conversely, when two variables are strongly correlated, the MIC approaches 1.

[0071] Furthermore, in step (6), the BPNN machine learning model has a three-layer structure, including an input layer, a hidden layer and an output layer, and the Tan-sigmoid and nonlinear LM algorithms are set as the transfer function and training function, respectively.

[0072] Beneficial effects: Compared with the prior art, the present invention has the following significant advantages: The present invention determines the main control parameters of the equivalent permeability of fracture networks by using the maximum information index method, and constructs a prediction model of the equivalent permeability of three-dimensional fracture networks based on the BPNN machine learning model. It comprehensively improves the shortcomings of high computational cost of discrete fracture models and the lack of consideration of topological connectivity in equivalent permeable medium models. It can efficiently and reliably convert fracture networks into corresponding equivalent media, and provides a simulation scheme with high prediction accuracy and low computational cost that reflects topological connectivity for the numerical simulation of underground flow fields in fractured strata at the site scale. Attached Figure Description

[0073] Figure 1 This is a flowchart of a numerical simulation method for estimating equivalent seepage flow in a three-dimensional fracture network based on topological properties, provided in an embodiment of the present invention.

[0074] Figure 2 This is a schematic diagram of the structure of the BPNN machine learning model in an embodiment of the present invention;

[0075] Figure 3 This is a schematic diagram illustrating the correlation analysis results between fracture network parameters and equivalent permeability in an embodiment of the present invention;

[0076] Figure 4 (a) is a target discrete fracture network diagram in an embodiment of the present invention. Figure 4 (b) is a diagram of the estimated equivalent permeation medium in an embodiment of the present invention;

[0077] Figure 5(a) is a top view of the target discrete fracture network in the vertical direction from 0 to 15 meters in an embodiment of the present invention. Figure 5 (b) is a top view of the estimated equivalent permeable medium in the vertical direction of the embodiment of the present invention, ranging from 0 to 15 meters.

[0078] Figure 6 (a) is a top view of the underground flow field of the target discrete fracture network in an embodiment of the present invention, in the vertical direction from 0 to 15 meters. Figure 6 (b) is a top view of the estimated equivalent permeable medium underground flow field in the vertical direction of the range of 0 to 15 meters in an embodiment of the present invention;

[0079] Figure 7 This is a top view of the target discrete fracture network, the estimated equivalent permeable medium, and the streamlines of the underground flow field of the estimated equivalent permeable medium within a vertical range of 0 to 15 meters in an embodiment of the present invention. Detailed Implementation

[0080] The invention will now be further described with reference to the accompanying drawings.

[0081] like Figure 1 As shown, this embodiment of the invention provides a numerical simulation method for estimating the equivalent seepage flow in a three-dimensional fracture network based on topological properties, comprising the following steps:

[0082] (1) Collect basic data on fracture networks and hydrogeological data, and construct a target discrete three-dimensional fracture network based on the basic data on fracture networks. The basic data on fracture networks include the distribution characteristics of fracture orientation, fracture diameter, and fracture density. The basic data on hydrogeology includes structural geological maps and hydrogeological maps.

[0083] (2) Construct a conceptual model of underground flow field of fracture network based on target discrete three-dimensional fracture network and basic hydrogeological data.

[0084] (3) Based on geological conditions, hydrogeological conditions and target area scale, determine the coverage range of fracture network parameters, and then use the Monte Carlo method to establish several small three-dimensional fracture networks by randomly generating fracture orientation to form an initial training sample set.

[0085] The attitude parameters of the small three-dimensional fracture network are as follows:

[0086] The crack diameter follows a power-law distribution, expressed as:

[0087] f(l)=βl -α (1)

[0088] Where f(l) is the probability density function of the crack diameter distribution; β is the probability density function based on the minimum value of the crack diameter l. min and maximum value l maxThe obtained normalization exponent; α is the power law exponent, set to 3.0;

[0089] The crack diameter is expressed as:

[0090]

[0091] Where R is a random number;

[0092] The fracture orientation follows the Fisher distribution:

[0093]

[0094] Where θ is the dip angle; θ0 is the average dip angle; φ is the dip direction; φ0 is the average dip direction; κ is the dispersion parameter, which takes a value of 0; and sinh is the hyperbolic sine function.

[0095] The fracture aperture λ ranges from 0.1 mm to 1 mm, conforming to a uniform random distribution. The fracture aperture differs from the fracture size by more than two orders of magnitude, and the effect of the fracture aperture size variation on seepage is negligible. Therefore, this invention does not consider the variation of fracture pore size and assumes that the fractures in each fracture network have the same aperture.

[0096] Based on the above parameter range and value selection method, the Monte Carlo method is used to generate the initial training sample set.

[0097] (4) Use graph theory to map the sample fracture network into an equivalent two-dimensional undirected graph, and solve for the equivalent permeability of the sample fracture network;

[0098] (4-1) Use graph theory to map the sample gap network into an equivalent two-dimensional undirected graph.

[0099] Treating the cracks as nodes V and the intersections between cracks as edges E, the 3D crack network is transformed into a 2D undirected graph G(V, E) as follows:

[0100] The three-dimensional fracture network F consists of N fractures f, represented as:

[0101] F={f i}, i = 1, ..., N

[0102] Define a bijective mapping Φ:

[0103] Φ:f i →v i

[0104] If there are two cracks f i with f j intersect, Then there is an edge in E connecting the corresponding vertex:

[0105]

[0106] Among them, (v i v j )∈E represents vertex v i and v j The edge e between ij , side e ij weight w ij Set as vertex v i and v j Length of the line of intersection between them:

[0107] w ij =Lengthf i ∩f j

[0108] Similarly, considering the flow direction, we can regard the inlet plane x0 as the source node s, and the outlet plane x L Considering node t as the target node, when the crack intersects with the inlet / outlet boundary,

[0109]

[0110]

[0111] Among them, e si For an edge between a source node and a non-source / sink node, e it The edges between non-source / sink nodes and the target node are thus transformed into a two-dimensional undirected graph.

[0112] When the 3D fracture network is sparse, the 2D undirected graph will be divided into several parts. Only the fracture sub-network corresponding to the connected subgraph containing source nodes and target nodes can effectively connect the inlet plane and the outlet plane and induce seepage. Therefore, the connected subgraph without source nodes and target nodes can be ignored. That is, the 2D undirected graph G can be pruned into the union of all connected subgraphs containing source nodes and target nodes except for edge nodes, i.e., a subset G' of the 2D undirected graph G. G' has a unique corresponding sub-network F' of the 3D fracture network F.

[0113] (4-2) Solving for the equivalent permeability of the sample fracture network

[0114] Solve for the equivalent permeability of the fracture network flow direction on the equivalent subgraph G' of the subnetwork F', based on the following assumptions:

[0115] i. The rock matrix is ​​incompressible and impermeable;

[0116] ii. The fluid is incompressible;

[0117] iii. The two planes of the crack are parallel and smooth;

[0118] iv. Fluid flow follows Darcy's law and exhibits laminar flow behavior.

[0119] Suppose there are M nodes in a two-dimensional undirected graph G, where each node represents the intersection of fractures, and node i is adjacent to node j. According to the cube law, the total flow Q between node i and node j is... j Represented as:

[0120]

[0121] Where k is the equivalent permeability; g is the gravitational acceleration; μ is the hydrodynamic viscosity; I is the hydraulic gradient; λ ij b is the aperture between crack i and crack j; b is the crack width; H i For the upstream head, H j For downstream head, L ij Let be the Euclidean distance between node i and node j;

[0122] According to the law of conservation of mass, the flow rate at node i is expressed as:

[0123]

[0124] The regional head distribution of the sample fracture network is represented as follows:

[0125]

[0126] Where, ΔH ij The difference in water head between upstream and downstream; w ij Let w be the equivalent water conductivity between node i and node j, and also the weight of the edge connecting node i and node j. When there is no connection between node i and node j, w is... ij The value is 0;

[0127] From this, we can obtain the adjacency matrix A and the degree matrix D of the vertices of a two-dimensional undirected graph G, defined as follows:

[0128] The equivalent permeability k of the entire sample fracture network region is expressed as:

[0129]

[0130] Among them, L zone Q is the region length of the sample crack network. si It is the inflow, H s It's the inlet head, H t It's the water outlet head.

[0131] (5) The maximum information index method is used to mine the main control parameters (or main control parameters) of equivalent permeability from several fracture network parameters, construct a sample set with the main control parameters as input and equivalent permeability as output, and normalize the input and output variables.

[0132] (5-1) The parameters of a fracture network include topological and geometric parameters.

[0133] Several important topology parameters are listed below:

[0134] Degree centrality represents the number of fractures that intersect with the fracture.

[0135] Eigenvector centrality is the iterative exponent of degree centrality, representing the importance of the fracture in the entire fracture network;

[0136] Clustering coefficient represents the diversity of pathways in a fracture network;

[0137] Betweenness centrality, which represents the number of the fracture in the shortest path of the network, plays an important role in controlling seepage.

[0138] Maximum flow is a feasible flow with the maximum capacity in a network with capacity, and it is an important topological feature of seepage flow.

[0139] The above topological parameters only represent the properties of a single fracture. Therefore, this invention takes their global average value and variance to reflect the seepage capacity and heterogeneity of the fracture network.

[0140] The anisotropy index fA is expressed as:

[0141]

[0142] Where i is the unit vector of the flow direction and n is the normal vector of the crack.

[0143] Geometric parameters include the number of fractures, fracture aperture, and fracture area P per unit volume. 32 and unit cell size.

[0144] (5-2) Use the maximum information index method to mine the main control parameters of the above topological parameters and a small number of geometric parameters.

[0145] The principle of the maximum information index method is as follows:

[0146] Suppose the dataset has Γ data points. Plot a scatter plot of the two attributes in a two-dimensional space, dividing it into a grid of size m×n along the horizontal and vertical axes, where m×n < B and B = Γ. 0.6 The grid with coordinates (x0, y0) contains γ data points, with a frequency P(x, y) = γ / Γ. Based on this, the normalized mutual information maximum grid resolution is calculated as the metric of MIC, i.e.:

[0147]

[0148] Here, X and Y are two variables whose correlation needs to be calculated.

[0149] The MIC value ranges from |0 to 1|. When two variables are independent, the MIC approaches 0; conversely, when two variables are strongly correlated, the MIC approaches 1.

[0150] (6) Based on the BPNN machine learning model, a prediction model of the equivalent permeability of the three-dimensional fracture network is constructed. The target discrete fracture network in step (1) is meshed. The three-dimensional matrix composed of the master control parameters is used as the input, and the three-dimensional matrix composed of the equivalent permeability, i.e. the equivalent medium, is used as the output. The hyperparameters of the model are adjusted to achieve high-precision simulation of input and output. Finally, a numerical model for simulating the underground flow field of the fracture network is obtained.

[0151] like Figure 2 As shown, the BPNN machine learning model includes an input layer, hidden layers, and an output layer. Each layer contains multiple neurons, and neurons in adjacent layers are interconnected through weighted connections, where w and V are the connection coefficients. The input layer serves as the data inlet, receiving input x, which is then passed to the hidden layer to obtain data h. The hidden layer calculates the output y from the data h. After minimizing the loss through the training function, backpropagation adjusts the connection coefficients in the network until the loss reaches the desired level or the training reaches the expected number of iterations. The BPNN machine learning model established in this invention has three layers, with Tan-sigmoid and the nonlinear LM (Levenberg-Marquardt) algorithm set as the transfer function and training function, respectively.

[0152] The BPNN machine learning model was trained using MATLAB software, and all sample analysis, spatial discretization, and numerical simulation in this invention were performed using MATLAB software.

[0153] The numerical simulation method for estimating equivalent seepage in a three-dimensional fractured network based on topological properties provided in this invention uses a BPNN machine learning model to predict the equivalent permeable medium containing topological properties in the formation. This greatly reduces the computational load of solving the underground flow field model using a discrete fractured network model, and provides an efficient and reliable simulation scheme for simulating seepage in fractured media under actual site conditions, with broad application prospects.

[0154] Here is a specific example.

[0155] This example demonstrates a numerical simulation of the equivalent permeability of a three-dimensional fractured network under actual site conditions. It comprehensively considers the topological properties of the fractured network and mainly includes two parts: 1) the mining of the main controlling factors of the equivalent permeability of the fractured network based on the maximum information index method; 2) the equivalent permeability prediction model based on the BPNN machine learning model.

[0156] 1) Mining the main controlling factors of equivalent penetration rate in fracture networks based on the maximum information index method

[0157] This example selects four typical orientations to prepare the sample set, and the specific discrete fracture network parameters are shown in Table 1. In scenario 1, the fracture diameter distribution is highly discrete; in scenario 2, the fracture diameter distribution is relatively uniform; in scenario 3, there is a set of conjugate fractures; in scenario 4, the fracture orientation conforms to a uniform random distribution, to emphasize the universality of the present invention.

[0158] The element side length in the simulation region is equal to the maximum crack diameter l. max To ensure accuracy, the fracture diameter changes proportionally with the size of the unit cell. Different simulation scenarios are set up to solve for the equivalent permeability of the fracture network in the main flow direction, and a database is established. This example controls the number of fractures between 5 and 4000 and the unit cell size between 0.55 and 100 meters, thus obtaining a large amount of data. After removing outliers using the k-means method, a database with a capacity of 5044 is obtained.

[0159] Table 1. Parameters of Discrete Fragment Network in Sample Set

[0160]

[0161] Where Δφ0 is the range of the difference between the dip angle and the mean dip angle, and Δθ0 is the range of the difference between the dip direction and the mean dip direction.

[0162] The maximum information index method was used to conduct a correlation analysis on the fracture network parameters and equivalent permeability of the sample set. Specific analysis results are as follows: Figure 3 As shown, Figure 3 In this context, lnk represents the logarithm of the equivalent permeability (the equivalent permeability is typically small, reaching at least 1E-20, therefore its logarithm is used for analysis and prediction). Fracture networks with different orientations share the same principal control parameters. Changes in unit cell size have little impact on parameter correlation, and geometric parameters show good correlation under different scenarios, while the correlation of topological parameters varies considerably. Fracture aperture, P... 32 The clustering coefficient variance exhibits a strong correlation with the equivalent penetration rate and a weak correlation with other parameters, making it the key control parameter for the prediction model. This simplifies model operation and improves accuracy.

[0163] 2) Equivalent penetration rate prediction model based on BPNN machine learning model

[0164] This example uses a synthetic fracture network, with parameters shown in Table 2. The BPNN machine learning model used in this example has a maximum training iteration count of 1000, a learning rate of 0.01, and a minimum training error of 0.00001. The input layer has 3 nodes, and the hidden layer has 5 nodes. The sample set is a cube with a side length of 15 meters, containing 16 to 336 discrete fractures. Based on the fracture survey results, 1050 small 3D fracture networks are generated. The dataset is divided into training and testing sets in an 8:2 ratio. The training set is used to train the BPNN machine learning model, and the testing set is used to verify its performance. The target region is a cube with a side length of 75 meters, containing 1000 discrete fractures. Based on the maximum fracture diameter, the target region is discretized into 125 sub-regions of 5×5×5 along the x, y, and z axes.

[0165] Table 2 Example Model Fractal Network Parameters

[0166]

[0167] By adjusting the number of training set samples, and using R... 2 The optimal training set size for the prediction model was determined using RMSE and MAPE, and the results are shown in Table 3. In all cases, the model performed well, and the predicted results showed a good fit with the calculated values, indicating that the method of this invention has high accuracy. The prediction model performed best with a dataset size of 1050, and the fit was slightly improved compared to the prediction results in other scenarios.

[0168] Table 3. Predicted Model Reliability Indicators

[0169]

[0170] like Figures 4 to 6 As shown, (a) is the discrete fracture model, and (b) is the estimated equivalent medium model. Figure 4 As can be seen, the prediction model converts the three-dimensional discrete fracture model of the site into an equivalent medium model. Figure 5 The model described is a discrete fracture model within the Z-axis range of 0–15 m, containing a total of 187 discrete fractures. The equivalent permeability is consistent with the distribution trend of the discrete fractures. At location (1,2), the fracture density of the mesh is high, exhibiting a high equivalent permeability. Conversely, the meshes at locations (1,3), (1,4), (5,2), and (5,3) show lower fracture density and lower equivalent permeability. Figure 6The steady-state seepage fields of the discrete fracture network model and the estimated equivalent medium model are presented respectively. Although the spatial discretization of the equivalent medium masks the strong inhomogeneity of the fractures to some extent, it effectively captures the flow obstruction and conduction phenomena caused by the uncertain fracture distribution at y=15m and y=45m, and the pressure front exhibits a distinct concave-convex shape. Figure 7 The flow field within the equivalent medium is shown. By comparing the discrete fractures and the equivalent medium, it is demonstrated that the streamline trend matches the fracture distribution. The fracture distribution is highly uneven near the inlet, with streamlines concentrated around the grid position (1,2). In the middle of the seepage path, the fracture distribution is more uniform, and the streamlines are more dispersed. The unevenness of the fracture distribution increases at the outlet, causing streamlines to move towards the region with denser fractures. This demonstrates the important role of fracture distribution in the seepage path and further proves the good ability of the simulation method of this invention to recreate the topological connectivity of the fracture network.

[0171] Overall, the trained BPNN machine learning model can effectively predict the equivalent medium distribution of a three-dimensional fracture network, realistically reflecting the inhomogeneity and topological connectivity of rock fractures, and achieving high prediction accuracy for the spatial distribution of underground seepage. Therefore, this prediction model can be used for three-dimensional visualization of the underground flow field of a three-dimensional fracture network under actual site conditions.

Claims

1. A numerical simulation method for estimating equivalent seepage flow in a three-dimensional fracture network based on topological properties, characterized in that, include: (1) Collect basic data on fracture networks and basic data on hydrogeology, and construct a target discrete three-dimensional fracture network based on the basic data on fracture networks; (2) Construct a conceptual model of the underground flow field of the fracture network based on the target discrete three-dimensional fracture network and basic hydrogeological data; (3) Based on geological conditions, hydrogeological conditions and target area scale, determine the coverage range of fracture network parameters, and then use the Monte Carlo method to establish several small three-dimensional fracture networks by randomly generating fracture orientation to form an initial training sample set. (4) Use graph theory to map the sample fracture network into an equivalent two-dimensional undirected graph, and solve for the equivalent permeability of the sample fracture network; (5) The maximum information index method is used to mine the main control parameters of equivalent permeability in several fracture network parameters, construct a sample set with the main control parameters as input and equivalent permeability as output, and normalize the input and output variables. (6) Based on the BPNN machine learning model, a prediction model of the equivalent permeability of the three-dimensional fracture network is constructed. The target discrete fracture network is gridded, and the three-dimensional matrix composed of the master control parameters is used as the input. The three-dimensional matrix composed of the equivalent permeability, i.e. the equivalent medium, is used as the output. The hyperparameters of the model are adjusted to achieve high-precision simulation of input and output. Finally, a numerical model for simulating the underground flow field of the fracture network is obtained.

2. The numerical simulation method according to claim 1, characterized in that, In step (1), the basic data of the fracture network includes the characteristics of fracture orientation distribution, fracture diameter distribution, and fracture density distribution.

3. The numerical simulation method according to claim 2, characterized in that, Basic hydrogeological data includes structural geological maps and hydrogeological maps.

4. The numerical simulation method according to claim 1, characterized in that, In step (3), the orientation parameters of the small three-dimensional fracture network are as follows: Crack diameter It follows a power-law distribution, represented as: in, Let be the probability density function of the crack diameter distribution; Based on the minimum value of the crack diameter and maximum value The obtained normalized index; The power-law exponent; Crack diameter Represented as: in, It is a random number; The fracture orientation follows the Fisher distribution: in, It is the angle of inclination; The average inclination angle; As a tendency; The tendency is towards the average; This is a parameter representing the degree of dispersion. It is a hyperbolic sine function; Crack aperture The values ​​range from 0.1 mm to 1 mm, conforming to a uniform random distribution; each fracture in the fracture network is set to have the same aperture.

5. The numerical simulation method according to claim 1, characterized in that, In step (4), graph theory is used to map the sample crack network into an equivalent two-dimensional undirected graph, including: Treating cracks as nodes The intersection line between the cracks is considered as an edge. Convert the three-dimensional fracture network into a two-dimensional undirected graph. The conversion method is as follows: Three-dimensional fracture network Depend on Strip crack Composition, represented as: , Define a bijective mapping : If two cracks and intersect, ,but There is an edge connecting the corresponding vertices: in, Represents vertices and The edge between ,side weight Set as vertex and Length of the line of intersection between them: Similarly, considering the flow direction, the inlet plane... Consider as source node , export plane Considered as the target node When the crack intersects with the inlet / outlet boundary: Among them, e si For an edge between a source node and a non-source / sink node, e it The edges between non-source / sink nodes and the target node are defined; thus, the 3D fracture network is transformed into a 2D undirected graph. When the 3D fracture network is sparse, the 2D undirected graph will be divided into several parts. Only the fracture sub-network corresponding to the connected subgraph containing source nodes and target nodes can effectively connect the inlet plane and the outlet plane, inducing seepage. Therefore, the connected subgraph without source nodes and target nodes is not considered; that is, the 2D undirected graph... The union of all connected subgraphs containing both the source and target nodes, excluding edge nodes, i.e., a two-dimensional undirected graph. subset of , There is a unique corresponding three-dimensional fracture network subnetwork .

6. The numerical simulation method according to claim 5, characterized in that, Solving for the equivalent permeability of the sample fracture network includes: in the subnetwork Equivalent subgraph The equivalent permeability in the flow direction of the fractured network is calculated based on the following assumptions: i. The rock matrix is ​​incompressible and impermeable; ii. The fluid is incompressible; iii. The two planes of the crack are parallel and smooth; iv. Fluid flow follows Darcy's law and exhibits laminar flow behavior.

7. The numerical simulation method according to claim 6, characterized in that, Assume a two-dimensional undirected graph There is There are nodes, where each node represents the intersection between fractures. With nodes Adjacent, according to the law of cubes, nodes With nodes Total flow between Represented as: in, Equivalent penetration rate; It is the acceleration due to gravity; The dynamic viscosity coefficient of the water flow; For hydraulic gradient; For cracks and cracks The degree of openness between them; The width of the crack; For the upstream water head, For downstream water head, For nodes With nodes The Euclidean distance between them; According to the law of conservation of mass, nodes The flow rate at the location is represented as: The regional head distribution of the sample fracture network is represented as follows: in, The difference in water head between upstream and downstream; For nodes With nodes The equivalent hydraulic conductivity between nodes is also the node With nodes The weight of the edges connecting the nodes, when the nodes With nodes When there is no connection between them, The value is 0; This yields a two-dimensional undirected graph. adjacency matrix The matrix of degree of vertices Defined as ; Equivalent permeability of the entire sample fracture network region Represented as: Among them, L zone It is the region length of the sample crack network. It is the inflow. It's the water inlet head. It's the water outlet head.

8. The numerical simulation method according to claim 1, characterized in that, In step (5), the fracture network parameters include topological and geometric parameters. The topological parameters include: Degree centrality represents the number of fractures that intersect with the fracture. Eigenvector centrality is the iterative exponent of degree centrality, representing the importance of the fracture in the entire fracture network; Clustering coefficient represents the diversity of pathways in a fracture network; Betweenness centrality, which represents the number of the fracture in the shortest path of the network, plays an important role in controlling seepage. Maximum flow is a feasible flow with the maximum capacity in a network with capacity, and it is an important topological feature of seepage flow. Anisotropy Index The expression is: in, The number of cracks. For crack aperture, The unit vector in the direction of flow. The crack normal vector; Geometric parameters include the number of fractures, fracture aperture, and fracture area P per unit volume. 32 and unit cell size.

9. The numerical simulation method according to claim 8, characterized in that, In step (5), the maximum information index method is used to mine the master control parameters, including: Assume the dataset has Given 10 data points, a scatter plot of the two attributes is drawn in two-dimensional space, divided along the horizontal and vertical axes. Dimensions of the grid , ; coordinates are The grid contains Data points, frequency Based on this, the maximum grid resolution of the normalized mutual information is calculated as... metric value: Where X and Y are the two variables whose correlation needs to be calculated; The range of values ​​is within When two variables are independent, The value approaches 0; conversely, when two variables are strongly correlated, Approaching 1.

10. The numerical simulation method according to claim 1, characterized in that, In step (6), the BPNN machine learning model has a three-layer structure, including an input layer, a hidden layer and an output layer. The Tan-sigmoid algorithm and the nonlinear LM algorithm are set as the transfer function and the training function, respectively.

Images