Evaluation Method of Hydraulic Fracturing Effect Based on FMM Accelerated Iterative Algorithm

Abstract

This invention discloses a method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm, comprising: constructing a three-dimensional fracture mesh model based on fracture geometric parameters and rock mechanics parameters; constructing a calculation equation for calculating the fracture width based on the displacement discontinuity method, wherein the calculation equation includes a coefficient matrix characterizing the elastic interaction between elements of the three-dimensional fracture mesh model; performing hierarchical compression of the coefficient matrix using a kernel-independent fast multi-level sub-method; obtaining rock mechanics parameters, initializing the mechanics parameters and the fracture state of all elements in the three-dimensional fracture mesh model; performing a double iteration to solve the fracture width; determining that the double iteration process has converged when the fracture state no longer changes and the relative change of the normal displacement vector is less than a preset convergence tolerance, and outputting the fracture width; and evaluating the fracture conductivity based on the fracture width.

Description

Technical Field

[0001] This invention relates to the field of oil and gas field development technology, and in particular to a method for evaluating the effect of hydraulic fracturing based on the FMM accelerated iterative algorithm. Background Technology

[0002] With the deepening development of unconventional oil and gas resources, accurate simulation of fracture closure behavior during hydraulic fracturing has become a crucial step in evaluating fracturing effectiveness and fracture conductivity. Fracture closure calculations involve solving large-scale linear systems, and traditional numerical simulation methods based on displacement discontinuity methods face challenges in computational efficiency and memory requirements when dealing with complex fracture networks. Particularly in three-dimensional fracture models, the storage and solution time complexity of the coefficient matrix increases quadratically with the number of mesh elements, limiting the widespread application of this method in engineering practice.

[0003] In existing technologies, to improve computational efficiency, a direct solver combined with matrix sparsification approximation is commonly used. However, while ensuring computational accuracy, these methods struggle to balance computational scale and often require simplified assumptions about crack geometry, leading to discrepancies between simulation results and actual conditions. Furthermore, traditional iterative solution methods are prone to convergence issues during crack state identification, particularly under complex conditions where crack opening and closing states alternate, making it difficult to guarantee algorithm stability and computational accuracy.

[0004] The Fast Multipole Method (FMM), as an efficient numerical algorithm, reduces computational complexity from O(N) to O(N) through hierarchical matrix compression techniques. 2 The computational cost has been reduced to O(NlogN), providing a new approach to solving large-scale boundary element problems. However, the application of traditional FMM in fracturing fracture simulation still faces a series of technical challenges, such as surrogate function selection, singularity handling, and boundary condition application. In particular, a mature and effective numerical implementation scheme has not yet been developed for the coupled solution of dynamic state changes during fracture closure.

[0005] Therefore, there is a lack of existing technologies that can balance computational efficiency, numerical accuracy, and physical rationality in fracture closure calculation. There is an urgent need to develop an efficient numerical simulation technology based on FMM acceleration and dual iterative algorithms to solve the key technical bottlenecks in large-scale fracture closure calculation and facilitate accurate evaluation of fracturing effects.

[0006] The above content is only used to help understand the technical solution of the present invention and does not represent an admission that the above content is prior art. Summary of the Invention

[0007] The main objective of this invention is to provide a method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm, which aims to solve or partially solve the above-mentioned problems.

[0008] To achieve the above objectives, this invention provides a method for evaluating the effectiveness of hydraulic fracturing based on the FMM accelerated iterative algorithm, comprising:

[0009] A three-dimensional crack mesh model is constructed based on the crack geometry parameters and rock mechanics parameters;

[0010] Based on the displacement discontinuity method, a computational equation for calculating crack width is constructed, which includes a coefficient matrix characterizing the elastic interaction between elements in the three-dimensional crack mesh model; and a kernel-independent fast multilevel sub-method is used to hierarchically compress the coefficient matrix.

[0011] Obtain rock mechanical parameters, initialize mechanical parameters and the crack state of all elements in the three-dimensional crack mesh model;

[0012] The crack width is solved by performing a double iteration on the computational equation. The double iteration includes an outer iteration and an inner iteration. The crack state is updated during the outer iteration, and the constrained linear system is solved during the inner iteration.

[0013] When the crack state no longer changes and the relative change of the normal displacement vector is less than the preset convergence tolerance, the double iteration process is considered to have converged, and the crack width is output.

[0014] Assess the flow conductivity of the crack based on its width.

[0015] Preferably, in the evaluation method for hydraulic fracturing effect based on the FMM accelerated iterative algorithm, the hierarchical compression of the coefficient matrix using a kernel-independent fast multi-level sub-method includes:

[0016] For near-field interactions in a 3D crack mesh model, the kernel function K ij As an element A of the coefficient matrix ij The corresponding calculation equation is: ;

[0017] For far-field interactions in a 3D crack mesh model, a kernel-independent fast multilevel sub-method is used to provide implicit matrix-vector multiplication operators. These multiplication operators are then used as computational equations. The multiplication operator is: ;

[0018] Where A is the coefficient matrix. ;

[0019] D n Let be the normal displacement vector to be solved, i.e., the crack width;

[0020] p is the fluid pressure vector inside the crack;

[0021] This represents the initial geostress vector;

[0022] F() is an implicit matrix-vector multiplication operator;

[0023] i represents the source element of the 3D crack mesh model, and j represents the field element of the 3D crack mesh model.

[0024] Preferably, in the evaluation method of hydraulic fracturing effect based on the FMM accelerated iterative algorithm, the interaction between source element i and field element j in the three-dimensional fracture mesh model is considered, and the kernel function K ij The calculation formula is:

[0025] ;

[0026] Among them, Cr=G / (4π(1-ν));

[0027] G is the shear modulus;

[0028] ν is Poisson's ratio;

[0029] b is the unit half-dimensional dimension, a = b = △x / 2;

[0030] (△x, △y, △z) is the relative position vector between the center points of the two units;

[0031] ;

[0032] ;

[0033] J 11 With J 21 This represents the basic integral kernel function defined under this mechanical model.

[0034] Preferably, in the evaluation method for hydraulic fracturing effect based on the FMM accelerated iterative algorithm, the step of providing implicit matrix-vector multiplication operators using a kernel-independent fast multilevel sub-method for far-field interactions in a three-dimensional fracture mesh model, and using the multiplication operators as computational equations, includes the following method for constructing the surrogate point set of the kernel-independent fast multilevel sub-method:

[0035] A basic proxy point set P0 is constructed by subdividing the spherical triangulation. A subset of 12 points is randomly selected from P0 as the basic direction vector.

[0036] Define the radius parameter geometry, R jFor each ∈ [1.5, 2.5], j = 1, 2, ..., k, a multi-radius spherical shell surrogate point set P is constructed by scaling each fundamental direction vector according to its radius value. multi ;

[0037] For the multi-radius spherical shell proxy point set P multi A 3D mesh alignment transformation is applied, and all transformed points form an aligned point set P. aligned This ensures the geometric consistency between the proxy points and the discrete mesh;

[0038] Aligned point set P aligned After deduplication, the final proxy point set P is formed. final .

[0039] Preferably, in the evaluation method of hydraulic fracturing effect based on FMM accelerated iterative algorithm, the rock mechanical parameters include rock elastic parameters, including elastic modulus E, Poisson's ratio v, and shear modulus in the step of obtaining rock mechanical parameters, initializing mechanical parameters and fracture state of all elements in the three-dimensional fracture mesh model.

[0040] Mechanical parameters include fluid pressure distribution, and the initialization formula is as follows:

[0041] ;

[0042] Where, p center The fluid pressure in the central region;

[0043] i is the source element of the three-dimensional crack mesh model;

[0044] p edge The fluid pressure at the edge / boundary;

[0045] For random fluctuations, ;

[0046] L center This defines the lateral boundary of the crack's central region.

[0047] In the 3D crack mesh model, the crack state (state) of all elements is a binary state vector, where state ∈ {0, 1}. N Initially, the crack state of all elements is set to open, and state(i) = 1.

[0048] Preferably, in the evaluation method for hydraulic fracturing effect based on the FMM accelerated iterative algorithm, the step of performing a double iteration to solve the fracture width in the calculation equation includes:

[0049] Iterate over each element of the 3D crack mesh model, if the source element i∈I open satisfy Then the crack state of source unit i will be changed from open to closed;

[0050] If source unit i∈I closed satisfy Then the crack state of source unit i will be changed from closed to open;

[0051] Among them, D n (i) represents the normal displacement of source element i;

[0052] ε d This is a preset displacement threshold;

[0053] The contact stress of source element i;

[0054] p(i) represents the fluid pressure;

[0055] This represents the initial ground stress;

[0056] ε s The preset stress threshold;

[0057] When source element i has a closed element, D n (i)=0, and the displacement constraints of the closed element are handled by the constraint projection function, and the modified linear system is solved by the preconditional conjugate gradient method.

[0058] Continue until the convergence condition is met.

[0059] Preferably, in the evaluation method for hydraulic fracturing effect based on the FMM accelerated iterative algorithm, the convergence condition in the step of "until convergence condition is met" simultaneously satisfies the following two conditions:

[0060] The crack state change converges, ||state （k） -state （k-1） ||=0;

[0061] The relative change in the normal displacement of the field element converges, as shown in the following formula:

[0062] ;

[0063] Where k is the number of iterations;

[0064] state （k） The crack state is obtained in the k-th iteration;

[0065] D n (k) This is the normal displacement obtained in the k-th iteration;

[0066] tol is the preset convergence tolerance.

[0067] Preferably, in the evaluation method for hydraulic fracturing effect based on the FMM accelerated iterative algorithm, the constraint projection function is defined as follows:

[0068] Diagonal projection matrix Its diagonal element S ii S ij From the state vector Sure:

[0069] ;

[0070] Constraining the projection function as an operator :R N →R N Its function is to process any input vector Perform the following operations:

[0071] ;

[0072] Where F() is the implicit matrix-vector multiplication operator;

[0073] state(i)=1 indicates that unit i is in an open state, and state(i)=0 indicates that unit i is in a closed state.

[0074] The present invention has at least the following beneficial effects:

[0075] The hydraulic fracturing effect evaluation method based on the FMM accelerated iterative algorithm provided by this invention constructs a three-dimensional fracture mesh model based on fracture geometric parameters and rock mechanics parameters; constructs a calculation equation for calculating the fracture width based on the displacement discontinuity method, wherein the calculation equation includes a coefficient matrix characterizing the elastic interaction between elements of the three-dimensional fracture mesh model; and performs hierarchical compression of the coefficient matrix using a kernel-independent fast multi-level sub-method; obtains rock mechanics parameters, initializes the mechanical parameters and the fracture state of all elements in the three-dimensional fracture mesh model; performs a double iteration to solve the fracture width, wherein the double iteration includes an outer iteration and an inner iteration. The outer iteration updates the fracture state, and the inner iteration solves the constrained linear system; when the fracture state no longer changes and the relative change of the normal displacement vector is less than the preset convergence tolerance, the double iteration process is considered to have converged, and the fracture width is output; the fracture conductivity is evaluated based on the fracture width. This can solve the key technical bottleneck in large-scale fracture closure calculation and provide a basis for fracturing effect evaluation.

[0076] The FMM acceleration technology used in this invention reduces the computational complexity from O(N) to O(N). 2 The time complexity is reduced to O(N log N), and the memory requirement is reduced from O(N log N). 2This invention reduces the computational complexity from O(N) to O(N), overcoming the computational bottleneck of traditional methods when dealing with large-scale crack meshes. It supports computation of meshes with tens of thousands of elements and provides technical support for simulating complex crack networks. The invention separates the precise calculation (near field) and approximate calculation (far field) of interactions between crack elements using a Fast Multipole (FMM) hierarchical algorithm. It also utilizes optimization techniques such as multi-radius spherical sampling to perform low-rank compression of far-field interactions, thereby reducing the computational complexity from O(N) of traditional methods. 2 The value is reduced to O(N log N).

[0077] Furthermore, the dual-iteration algorithm designed in this invention solves the convergence problem in crack closure calculation. The strategy of separating the outer state identification and the inner equation solution ensures physical rationality and numerical stability, and is applicable to complex working conditions with changing crack states.

[0078] Furthermore, the multi-radius spherical surrogate function and mesh alignment technique proposed in this invention avoid the singularity problem of traditional planar surrogate functions, improve the numerical accuracy of far-field approximation, and reduce the amount of computation by utilizing symmetry.

[0079] Furthermore, compared to traditional fracture closure analysis methods, this invention is based on the complete displacement discontinuity method theory, which can simulate fracture closure behavior under complex geometries and is applicable to the evaluation of fracturing effects in various unconventional reservoirs, thus having stronger applicability.

[0080] Furthermore, this invention introduces projection techniques and PCG iteration into the solution of constrained linear systems, effectively handling the displacement constraints of closed elements and ensuring the consistency between numerical solutions and physical constraints, which forms the basis for evaluating the conductivity of cracks. Attached Figure Description

[0081] Figure 1 A schematic diagram illustrating the evaluation method for hydraulic fracturing effect based on the FMM accelerated iterative algorithm provided by this invention;

[0082] Figure 2 A comparison chart of computation time between the direct method and the FMM method;

[0083] Figure 3 A comparison chart of memory requirements between the direct method and the FMM method;

[0084] Figure 4 The graph shows the relationship between the speedup ratio and the number of cells for the direct method and the FMM method.

[0085] Figure 5 This is a comparison chart of the seam widths of the direct method and the FMM method.

[0086] The realization of the objective, functional features and advantages of the present invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0087] In this embodiment of the invention, the term "and / or" describes the relationship between associated objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, or B existing alone. The character " / " generally indicates that the preceding and following associated objects have an "or" relationship.

[0088] It should be noted that the terms "first," "second," etc., in the specification, claims, and drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence.

[0089] In this embodiment of the invention, the term "multiple" refers to two or more, and other quantifiers are similar.

[0090] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the various embodiments of the present invention will be described in detail below with reference to the accompanying drawings. However, those skilled in the art will understand that many technical details are presented in the embodiments of the present invention to facilitate a better understanding of the invention. However, the technical solutions claimed in the present invention can be implemented even without these technical details and various variations and modifications based on the following embodiments. The division of the following embodiments is for ease of description and should not constitute any limitation on the specific implementation of the present invention. The various embodiments can be combined with and referenced by each other without contradiction.

[0091] This invention provides a method for evaluating the effectiveness of hydraulic fracturing based on the FMM accelerated iterative algorithm, such as... Figure 1 As shown, Figure 1 The diagram illustrates the process of evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm provided by this invention.

[0092] In step S100, a three-dimensional fracture mesh model is constructed based on the fracture geometry parameters and rock mechanics parameters. The fracture geometry parameters include fracture half-length, fracture half-height, and spatial step size. The rock mechanics parameters include rock elastic modulus and Poisson's ratio. Based on the fracture geometry parameters and rock mechanics parameters, the fracture surface is discretized in three-dimensional space to generate a planar mesh model containing N elements.

[0093] It should be noted that the spatial step size controls the calculation accuracy. The smaller the step size, the denser the grid, and the higher the calculation accuracy, but the amount of computation also increases accordingly.

[0094] In step S200, a calculation equation for calculating the crack width is constructed based on the displacement discontinuity method. The calculation equation includes a coefficient matrix that characterizes the elastic interaction between elements in the three-dimensional crack mesh model. The coefficient matrix is ​​then hierarchically compressed using a kernel-independent fast multilevel sub-method.

[0095] The governing equations for the Discontinuous Displacement Method (DDM) are: .

[0096] A is the coefficient matrix. ; Elements A of the coefficient matrix ij From kernel function K ij The calculated values ​​characterize the elastic interactions between the elements. K ij The normal stress induced at the center point of source element i by a unit normal displacement discontinuity on source element j in an infinitely large elastic medium is defined, which is the theoretical basis of the displacement discontinuity method.

[0097] D n ∈R N Let be the normal displacement vector to be solved, i.e., the crack width;

[0098] p∈R N This is the fluid pressure vector within the crack;

[0099] ∈R N This is the initial geostress vector.

[0100] Considering the interaction between source element i and field element j in a 3D crack mesh model, the kernel function K ij The calculation formula is:

[0101] ;

[0102] Among them, Cr=G / (4π(1-ν));

[0103] G is the shear modulus;

[0104] ν is Poisson's ratio;

[0105] b is the unit half-dimensional dimension, a = b = △x / 2;

[0106] (△x, △y, △z) is the relative position vector between the center points of the two units;

[0107] ;

[0108] ;

[0109] J 11 With J 21 This represents the basic integral kernel function defined under this mechanical model.

[0110] Directly processing dense matrices A typically presents computational bottlenecks. This invention employs a kernel-independent fast multilevel sub-method to accelerate this process.

[0111] Specifically, for near-field interactions in a 3D crack mesh model, the kernel function K... ij As an element A of the coefficient matrix ij The corresponding calculation equation is: ;

[0112] For far-field interactions in a 3D crack mesh model, a kernel-independent fast multilevel sub-method is used to provide implicit matrix-vector multiplication operators. These multiplication operators are then used as computational equations. The multiplication operator is: ;

[0113] Where A is the coefficient matrix. ;

[0114] D n Let be the normal displacement vector to be solved, i.e., the crack width;

[0115] p is the fluid pressure vector inside the crack;

[0116] This represents the initial geostress vector;

[0117] F() is an implicit matrix-vector multiplication operator;

[0118] i represents the source element of the 3D crack mesh model, and j represents the field element of the 3D crack mesh model.

[0119] More specifically, an octree spatial partitioning structure covering all crack units is constructed. The number of units contained in a leaf node is controlled by the parameter occ, where occ is set to 256. Due to the merging and sorting of row and column points, the actual number of units contained in each leaf node is approximately 128. For near-field interactions (direct neighbor nodes in the tree structure), the matrix elements of the coefficient matrix A are obtained by directly calculating the kernel function. For far-field interactions, low-rank compression is performed using interpolation decomposition: the interaction matrix is ​​calculated using the surrogate point set, and interpolation decomposition is performed to obtain the skeleton points and redundant points, thus constructing a hierarchical compression structure. The output of operator F() is obtained through the processes of upward scanning (compression), interaction calculation, and downward scanning (refinement).

[0120] It should be noted that when occ is 256, the FMM algorithm achieves optimal performance in the fracture closure problem of hydraulic fracturing.

[0121] The method for constructing the proxy point set of the kernel-independent fast multi-level sub-method includes steps S210 to S240.

[0122] In step S210, a basic proxy point set P0 is constructed by subdividing the spherical triangulation. A subset of 12 points is randomly selected from P0 as the basic direction vector. ,in, Let be a point on the unit sphere. Through random sampling... Select a subset containing 12 points , which serves as the fundamental direction vector.

[0123] Define the radius parameter geometry at step S220, R j For each ∈ [1.5, 2.5], j = 1, 2, ..., k, a multi-radius spherical shell surrogate point set P is constructed by scaling each fundamental direction vector according to its radius value. multi .

[0124] Define a set of radius parameters {R1,R2,……,R k}, where R j ∈[1.5,2.5], j=1,2,……,k. Twelve radius values ​​are selected at equal intervals, and this basic proxy point set is used... That is, each basic direction vector is scaled according to its radius value to form 144 surrogate points (12 basic directions × 12 radii), thus constructing a multi-radius spherical shell surrogate point set P. multi .

[0125] .

[0126] At step S230, the multi-radius spherical shell proxy point set P multi A 3D mesh alignment transformation is applied, and all transformed points form an aligned point set P. aligned To ensure geometric consistency between the surrogate points and the discrete mesh, let Δh be the spatial discrete step size of the crack plane. For the multi-radius spherical shell surrogate point set P... multi For each point in the equation, q = (q x ,q y ,q z )∈P multi Perform a grid alignment transformation as follows:

[0127] ;

[0128] ;

[0129] .

[0130] Here, `round()` is the rounding function. All transformed points form the aligned point set P. aligned .

[0131] The surrogate point technique of general FMM is designed for smooth field problems such as electrostatic fields. However, the DDM kernel function involved in hydraulic fracturing has nonlinearity, high gradient and singularity in the near field, which cannot be accurately approximated by general methods.

[0132] This invention employs a targeted multi-radius spherical sampling technique: this technology improves the approximate accuracy of the FMM for the DDM kernel function in the near-far field transition region. The DDM kernel function exhibits nonlinearity and high gradient changes in the near field, making it difficult to accurately capture using traditional single spherical surrogate points. Multi-radius spherical sampling, by covering a spatial range, avoids singularities and provides a high-precision interpolation basis. Combining "multi-radius spherical sampling" with "mesh alignment technology" to optimize FMM performance is a solution to the unique challenge of near-field accuracy and stability of the DDM kernel function.

[0133] The point set P after alignment at step S240 aligned After deduplication, the final proxy point set P is formed. final .

[0134] From the aligned point set P aligned Remove all duplicate coordinate points to form the final FMM surrogate point set P. final This is used for subsequent fast multipole calculations.

[0135] Multi-radius spherical sampling is introduced into the DDM simulation of hydraulic fracturing fractures. By constructing a spherical shell instead of a sphere, it systematically improves the approximate accuracy of the kernel function throughout the far field, which is fundamental to ensuring the accuracy of fracture closure determination. Mesh alignment technology forces the mathematical surrogate point set to be geometrically aligned with the physically discrete mesh, enhancing the numerical stability and computational efficiency of the algorithm, and is key to achieving efficient O(Nlog N) computation.

[0136] In step S300, rock mechanical parameters are obtained, and the mechanical parameters and the crack state of all elements in the three-dimensional crack mesh model are initialized. Rock mechanical parameters include rock elastic parameters, including elastic modulus E, Poisson's ratio v, and shear modulus.

[0137] Mechanical parameters include fluid pressure distribution, and the initialization formula is as follows:

[0138] ;

[0139] Where, p center The fluid pressure in the central region;

[0140] i is the source element of the three-dimensional crack mesh model;

[0141] p edge The fluid pressure at the edge / boundary;

[0142] For random fluctuations, ;

[0143] L center This defines the lateral boundary of the crack's central region.

[0144] In the 3D crack mesh model, the crack state (state) of all elements is a binary state vector, where state ∈ {0, 1}. N Initially, the crack state of all elements is set to open, and state(i) = 1.

[0145] In step S400, the crack width is solved by double iteration of the calculation equation. The double iteration includes outer iteration and inner iteration. The crack state is updated during the outer iteration, and the constrained linear system is solved during the inner iteration.

[0146] Specifically, step S400 includes steps S410 to S440.

[0147] In step S410, each element of the 3D crack mesh model is iterated. If the source element i∈I open satisfy Then the crack state of source element i will be changed from open to closed.

[0148] Let I open with I closed These represent the sets of cell indices that are in the open and closed states in the current iteration step, respectively. In this embodiment, we take... m.

[0149] If the source element i∈I at step S420 closed satisfy Then the crack state of source unit i will be changed from closed to open;

[0150] Among them, D n (i) represents the normal displacement of source element i;

[0151] ε d This is a preset displacement threshold;

[0152] The contact stress of source element i;

[0153] p(i) represents the fluid pressure;

[0154] This represents the initial ground stress;

[0155] ε s This is a preset stress threshold. In this embodiment, .

[0156] In step S430, when source element i has a closed element, D is... n (i)=0, and the displacement constraints of the closed element are handled by the constraint projection function. The modified linear system is solved using the preconditional conjugate gradient method. The constraint projection function is defined as follows:

[0157] Diagonal projection matrix Its diagonal element S ii S ij From the state vector Sure:

[0158] ;

[0159] Constraining the projection function as an operator :R N →R N Its function is to process any input vector Perform the following operations:

[0160] ;

[0161] Where F() is the implicit matrix-vector multiplication operator;

[0162] state(i)=1 indicates that unit i is in an open state, and state(i)=0 indicates that unit i is in a closed state.

[0163] Operator Ensure that the displacement of the closed element is forced to zero; remove the corresponding equation of the closed element; accelerate the calculation of far-field effects through FMM.

[0164] By fusing the projection matrix S with the FMM implicit operator F, displacement constraints can be handled without explicitly constructing, modifying, or storing the dense coefficient matrix A.

[0165] At step S440, the convergence condition is met until it is satisfied. The convergence condition, in this case, simultaneously satisfies the following two conditions:

[0166] The crack state change converges, ||state （k） -state （k-1） ||=0; the first norm of the state vector (the sum of the absolute values ​​of each component) is zero, indicating that the open / closed states of all elements no longer change.

[0167] The relative change in the normal displacement of the field element converges, as shown in the following formula:

[0168] ;

[0169] Where k is the number of iterations;

[0170] state （k） The crack state is obtained in the k-th iteration;

[0171] D n (k)This is the normal displacement obtained in the k-th iteration;

[0172] `tol` is the preset convergence tolerance. The relative change of the normal displacement vector is less than the preset convergence tolerance `tol`, and the maximum number of iterations is set to `k`. max =1000, where tol is 10 −8 .

[0173] This invention, through dual iteration, ensures convergence to the most accurate solution for the crack width.

[0174] In step S500, when the crack state no longer changes and the relative change of the normal displacement vector is less than the preset convergence tolerance, the convergence of the dual iteration process is determined, and the crack width is output.

[0175] In step S600, the flow conductivity of the fracture is assessed based on its width. Assessing flow conductivity based on fracture width is a standard method and will not be detailed here. A wider fracture generally indicates a greater flow conductivity.

[0176] Example:

[0177] Numerical verification was performed using a three-dimensional fracture model. The case parameters were set as follows: fracture half-length 150m, half-height 30m, rock elastic modulus 30.0 GPa, Poisson's ratio 0.2, and initial in-situ stress 10 MPa. The fluid pressure distribution within the fracture was non-uniform, with a pressure of 10 MPa in the central region (|x|≤100m) and 9.5 MPa in the edge region (|x|>100m). Random perturbations were superimposed to simulate actual working conditions. The spatial step size dx was set to 3m, 2m, 1m, and 0.5m, corresponding to 2121, 4681, 18361, and 72721 mesh elements, respectively. The maximum number of iterations was set to 1000, and the convergence tolerance was 10. −8 .

[0178] Solution: According to the evaluation method of hydraulic fracturing effect based on FMM accelerated iterative algorithm provided by this invention, firstly, a fracture geometric model is constructed and mechanical parameters are initialized. Then, the coefficient matrix is ​​constructed using the FMM accelerated displacement discontinuity method, where the surrogate function adopts multi-radius spherical sampling technology with a radius range of [1.5, 2.5] and the number of surrogate points p=144. The fracture closure problem is solved through a dual iterative algorithm, with the outer iteration updating the fracture state and the inner iteration using PCG to solve the constrained linear system.

[0179] Computation time, for example Figure 2As shown, as the number of mesh elements increases from 2121 to 72721, the computation time of the FMM method increases from 0.1s to 4.1s, while the direct method already has a long computation time when the number of elements exceeds 4681. At a mesh size of 18361 elements, the speedup of the FMM method compared to the direct method is (…). Figure 4 The efficiency is 45 times higher, demonstrating the advantages of this invention in large-scale computing.

[0180] Memory requirements, for example Figure 3 As shown, the memory consumption of the FMM method increases approximately linearly with the grid size, from 6.4MB at 2121 cells to 367.4MB at 72721 cells; while the memory consumption of the direct method increases quadratically, reaching 175.3MB at 4681 cells, and theoretically exceeding 2.5GB at 18361 cells, which exceeds the capacity of a conventional computer.

[0181] seam width for example Figure 5 As shown, to verify the computational accuracy of the FMM acceleration method, this invention compares the crack width distributions obtained by the direct DDM method and the FMM acceleration method. By extracting crack width data along the crack centerline (y=0) and conducting comparative analysis, it was found that the crack width distribution curves calculated by the two methods basically overlap, which fully demonstrates that the FMM acceleration method improves computational efficiency while maintaining computational accuracy.

[0182] This example verifies the efficiency and feasibility of the method of the present invention in large-scale crack closure calculation, and provides technical support for practical engineering applications.

[0183] Obviously, the embodiments described above are merely some, not all, embodiments of the present invention. Based on the embodiments of the present invention, those skilled in the art can make other variations or modifications without creative effort, and all such variations or modifications should fall within the scope of protection of the present invention.

Claims

1. A method for evaluating the effect of hydraulic fracturing based on the FMM accelerated iterative algorithm, characterized in that, The method comprises the following steps: According to the crack geometric parameters and rock mechanical parameters, a three-dimensional fracture grid model is constructed; Based on the displacement discontinuity method, a calculation equation for calculating the fracture width is constructed, wherein the calculation equation includes a coefficient matrix representing the elastic interaction between the units of the three-dimensional fracture grid model; and the coefficient matrix is compressed hierarchically using a kernel-independent fast multilevel cell method; Obtain rock mechanical parameters, initialize mechanical parameters and the crack state of all units in the three-dimensional fracture grid model; Solve the crack width by double iteration of the calculation equation, wherein the double iteration includes outer iteration and inner iteration, the crack state is updated during the outer iteration, and the constrained linear system is solved during the inner iteration; When the crack state no longer changes and the relative change of the normal displacement vector is less than the preset convergence tolerance, it is determined that the double iteration process converges, and the crack width is output; According to the crack width, the fracture conductivity is evaluated; The hierarchical compression of the coefficient matrix using the kernel-independent fast multilevel cell method comprises: For the near-field interaction in the three-dimensional fracture network model, the kernel function Kij is taken as the element Aij of the coefficient matrix, and the corresponding calculation equation is: ; For far-field interactions in three-dimensional fracture mesh models, a kernel-independent fast multilevel cell method is used to provide an implicit matrix-vector multiplication operator, which is used as the computational equation, where ; wherein A is a coefficient matrix, ; Dn is the normal displacement vector to be solved, i.e. the crack width; p is the fluid pressure vector in the crack; is the initial stress vector; F() is an implicit matrix-vector multiplication operator; i is a source unit of the three-dimensional fracture grid model, and j is a field unit of the three-dimensional fracture grid model.

2. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 1, wherein, Considering the interaction between the source unit i and the field unit j in the three-dimensional fracture grid model, the calculation formula of the kernel function Kij is: ； where Cr=G / (4π(1-ν)); G is the shear modulus; ν is the Poisson's ratio; b is the unit half size, a=b=△x / 2; (△x,△y,△z) is the relative position vector between the centers of two units; ； ； J11 and J21 represent the basic integral kernel function defined under the mechanical model.

3. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 1, wherein, For far-field interaction in the three-dimensional fracture grid model, the kernel-independent fast multilevel cell method provides an implicit matrix-vector multiplication operator, and the construction method of the proxy point set of the kernel-independent fast multilevel cell method includes: A basic proxy point set P0 is constructed by spherical triangle subdivision, and a subset of 12 points is randomly sampled from P0 as a basic direction vector; Define the radius parameter geometry Rj∈[1.5,2.5], j=1,2, …, k, and construct a multi-radius spherical shell proxy point set Pmulti by scaling each basic direction vector according to the radius value; A three-dimensional grid alignment transformation is performed on the multi-radius spherical shell proxy point set Pmulti, and all the transformed points form an aligned point set Paligned, ensuring the geometric consistency of the proxy points and the discrete grid; The final proxy point set Pfinal is formed by removing duplicates in the aligned point set Paligned.

4. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 1, wherein, In the step of obtaining rock mechanical parameters, initializing mechanical parameters and the crack state of all units in the three-dimensional fracture grid model, the rock mechanical parameters include rock elastic parameters, including elastic modulus E, Poisson's ratio v, and shear modulus; The mechanical parameters include fluid pressure distribution, and the initialization formula is as follows: ； where pcenter is the fluid pressure of the central region; i is a source unit of the three-dimensional fracture grid model; pedge is the fluid pressure of the edge / boundary part; for the random fluctuation term, ; Lcenter is the lateral range limit of the crack center region; The crack state state of all units in the three-dimensional fracture grid model is a binary state vector, state∈{0,1}N, and the crack state state of all units is set to open at the initial time, state(i)=1.

5. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 4, characterized in that, The step of solving the fracture width by double iteration of the calculation equation comprises: Each element of the three-dimensional fracture mesh model is iterated, if the source element i∈Iopen satisfies the fracture state state of the source element i is converted from opening to closing; If the source unit i e Iclosed satisfies then the crack state state of the source unit i is changed from closed to open; Wherein, Dn(i) is the normal displacement of the source unit i; εd is a preset displacement threshold; contact stress for source unit i; p(i) is the fluid pressure; is the initial stress; εs is a preset stress threshold; When the source unit i has a closed unit, Dn(i)=0, and the displacement constraint of the closed unit is processed by a constraint projection function, a preconditioned conjugate gradient method is used to solve the modified linear system; Until the convergence condition is met.

6. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 5, wherein, The convergence condition in the step of until the convergence condition is met meets the following two conditions simultaneously: The crack state state changes converges, ||state(k)-state(k-1)||=0; The relative change of the normal displacement of the field unit converges, and the formula is as follows: ； Wherein, k is the iteration number; state(k) is the crack state obtained by the kth iteration; Dn(k) is the normal displacement obtained by the kth iteration; tol is a preset convergence tolerance.

7. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 5, characterized in that, The definition of the constraint projection function is as follows: Diagonal projection matrix whose diagonal elements Sii, Sij are determined from the state vector S = [Sii, Sij]T The constraint projection function as an operator : RN→ RN, which acts on any input vector performs the following operation: ； Wherein, F() is an implicit matrix vector multiplication operator; state(i)=1 indicates that the unit i is in an open state, and state(i)=0 indicates that the unit i is in a closed state.

Images