Hydrofracture effect evaluation method based on FMM accelerated iterative algorithm

Abstract

The invention discloses a hydrofracture effect evaluation method based on an FMM accelerated iterative algorithm, and the method comprises the steps: constructing a three-dimensional fracture mesh model according to fracture geometric parameters and rock mechanical parameters; on the basis of a displacement discontinuity method, a calculation equation used for calculating the crack width is constructed, and the calculation equation comprises a coefficient matrix representing elastic interaction between units of the three-dimensional crack grid model; performing hierarchical compression on the coefficient matrix by using a kernel-independent fast multilevel sub-method; rock mechanical parameters, initialized mechanical parameters and fracture states of all units in the three-dimensional fracture grid model are obtained; performing double iteration on the calculation equation to solve the crack width; when the crack state does not change any more and the relative change of the normal displacement vector is smaller than a preset convergence tolerance, it is determined that the dual iteration process is converged, and the crack width is output; and evaluating the fracture conductivity according to the fracture width.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of oil and gas field development, and particularly relates to a method for evaluating hydraulic fracturing effect based on a fast multipole method (FMM) accelerated iterative algorithm. BACKGROUND

[0002] With the deepening of the development of unconventional oil and gas resources, the accurate simulation of the fracture closure behavior in the hydraulic fracturing process has become a key link for evaluating the fracturing effect and the fracture conductivity. The fracture closure calculation involves the solution of a large-scale linear system. The traditional numerical simulation method based on the displacement discontinuity method faces the challenges of computational efficiency and memory requirement when dealing with complex fracture networks. In particular, in a three-dimensional fracture model, the storage of the coefficient matrix and the time complexity of the solution increase in square with the number of grid cells, which restricts the wide application of the method in engineering practice.

[0003] In the prior art, in order to improve the calculation efficiency, a direct solver combined with matrix sparsification approximation method is usually used, but this kind of method is difficult to balance the calculation scale while ensuring the calculation accuracy, and often needs to make simplifying assumptions on the fracture geometry, resulting in deviations between the simulation results and the actual situation. On the other hand, the traditional iterative solution method is prone to convergence problems in the fracture state identification process, especially in the complex working conditions where the fracture opening and closing states change alternately, the stability and calculation accuracy of the algorithm are difficult to guarantee.

[0004] As a kind of efficient numerical algorithm, the fast multipole method (FMM) reduces the computational complexity from O(N 2 ) to O(NlogN) through hierarchical matrix compression technology, providing a new idea for solving large-scale boundary element problems. However, the application of traditional FMM in fracture simulation still faces a series of technical difficulties such as proxy function selection, singularity processing and boundary condition application. Especially in the coupled solution of dynamic state changes in the fracture closure process, there is still no mature and effective numerical implementation scheme.

[0005] Therefore, there is a lack of a fracture closure calculation method that can balance the calculation efficiency, numerical accuracy and physical reasonableness in the prior art, and it is urgent to develop an efficient numerical simulation technology based on FMM acceleration and double iterative algorithm to solve the key technical bottlenecks in large-scale fracture closure calculation, so as to accurately evaluate the fracturing effect.

[0006] The above content is only used to assist in understanding the technical solutions of the present application, and does not represent the acknowledgement of the above content as prior art. SUMMARY

[0007] The main purpose of the present application is to provide a hydraulic fracturing effect evaluation method based on FMM accelerated iterative algorithm, aiming to solve or partially solve the above problems.

[0008] To achieve the above-mentioned purpose, the present application provides a hydraulic fracturing effect evaluation method based on FMM accelerated iterative algorithm, comprising: According to the crack geometric parameters and rock mechanics parameters, a three-dimensional crack grid model is constructed; Based on the displacement discontinuity method, a calculation equation for calculating the crack width is constructed, wherein the calculation equation includes a coefficient matrix representing the elastic interaction between the units of the three-dimensional crack grid model; and the hierarchical compression of the coefficient matrix is performed by using the kernel-independent fast multi-level cell method; Obtain the rock mechanics parameters, initialize the mechanics parameters and the crack state of all units in the three-dimensional crack 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 crack conductivity is evaluated. Preferably, in the hydraulic fracturing effect evaluation method based on FMM accelerated iterative algorithm, the hierarchical compression of the coefficient matrix by using the kernel-independent fast multi-level cell method comprises: For the near-field interaction in the three-dimensional crack grid model, the kernel function K ij is taken as the element A ij of the coefficient matrix, ; For the far-field interaction in the three-dimensional crack grid model, the kernel-independent fast multi-level cell method is used to provide an implicit matrix vector multiplication operator, and the multiplication operator is taken as the calculation equation, and the multiplication operator is ; wherein A is the coefficient matrix, ; D n 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 ground stress vector; F() is the implicit matrix vector multiplication operator; i is the source unit of the three-dimensional crack grid model, and j is the field unit of the three-dimensional crack grid model.

[0009] Preferably, in the method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm, the interaction between the source unit i and the field unit j in the three-dimensional fracture grid model is considered, and the kernel function K ij The calculation formula of the kernel function K ; Wherein, Cr=G / (4π (1-v) ) ; G is the shear modulus; v is the Poisson's ratio; b is the half size of the unit, a=b=△x / 2; (△x,△y,△z) is the relative position vector between the centers of two units; ; ; J 11 And J 21 represent the basic integral kernel function defined under the mechanical model.

[0010] Preferably, in the method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm, for the far-field interaction in the three-dimensional fracture grid model, the kernel-independent fast multi-level cell method is used to provide an implicit matrix-vector multiplication operator, and the construction method of the proxy point set of the kernel-independent fast multi-level cell method in the step of calculating the equation includes: A basic proxy point set P0 is constructed by spherical triangle subdivision, and a subset containing 12 points is randomly sampled from P0 as a basic direction vector; A radius parameter geometry R j ∈[1.5,2.5], j=1,2, …, k, and a multi-radius spherical shell proxy point set P multi is constructed by scaling each basic direction vector according to the radius value. The multi-radius spherical shell proxy point set P multi is transformed by a three-dimensional grid alignment transformation, and all the transformed points form an aligned point set P aligned , which ensures the geometric consistency of the proxy points and the discrete grid; The aligned point set P aligned is de-duplicated to form a final proxy point set P final .

[0011] Preferably, in the method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm, in the step of obtaining the rock mechanics parameters, initializing the mechanics parameters and the fracture state of all units in the three-dimensional fracture grid model, the rock mechanics parameters include the rock elastic parameters, including the elastic modulus E, the Poisson's ratio v, and the shear modulus. The mechanical parameter includes fluid pressure distribution, and the initialization formula is as follows: ; Wherein, p center is the fluid pressure of the central region; i is the source element of the three-dimensional fracture grid model; p edge is the fluid pressure of the edge / boundary part; is a random fluctuation term, ; L center is the lateral range limit of the fracture central region; The fracture state state of all elements in the three-dimensional fracture grid model is a binary state vector, state∈{0,1} N , and the fracture state state of all elements is set to open at the initial time, state(i)=1.

[0012] Preferably, in the evaluation method of the hydraulic fracturing effect based on the FMM accelerated iterative algorithm, the step of solving the fracture width by double iteration of the calculation equation comprises: Iterate each element of the three-dimensional fracture grid model, if the source element i∈I open satisfies , the fracture state state of the source element i is converted from open to closed; If the source element i∈I closed satisfies , the fracture state state of the source element i is converted from closed to open; Wherein, D n (i) is the normal displacement of the source element i; ε d is a preset displacement threshold; is the contact stress of the source element i; p(i) is the fluid pressure; is the initial ground stress; ε s is a preset stress threshold; When the source element i has a closed element, D n (i)=0, and the displacement constraint of the closed element 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 satisfied.

[0013] Preferably, in the method for evaluating hydraulic fracturing effect based on the FMM accelerated iterative algorithm, the convergence condition is met simultaneously with the following two conditions: The crack state state 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 in the kth iteration; D n (k) is the normal displacement obtained in the kth iteration; tol is a preset convergence tolerance.

[0014] Preferably, in the method for evaluating hydraulic fracturing effect based on the FMM accelerated iterative algorithm, the definition of the constraint projection function is as follows: The diagonal projection matrix , the diagonal elements S ii , S ij are determined by the state vector : ; The constraint projection function is an operator : R N → R N , which performs the following operation on any input vector : ; 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.

[0015] The present application has at least the following beneficial effects: The application provides the evaluation method for the hydraulic fracturing effect based on the FMM accelerated iterative algorithm, a three-dimensional fracture grid model is constructed according to fracture geometric parameters and rock mechanics parameters; a calculation equation for calculating fracture width is constructed based on a displacement discontinuity method, wherein the calculation equation comprises a coefficient matrix representing elastic interaction between units of the three-dimensional fracture grid model; and the coefficient matrix is hierarchically compressed by using a kernel-independent fast multi-level cell method; rock mechanics parameters are acquired, and the mechanics parameters and the fracture state of all units in the three-dimensional fracture grid model are initialized; the calculation equation is solved by double iteration to obtain the fracture width, wherein the double iteration comprises outer iteration and inner iteration, the fracture state is updated during the outer iteration, and a constrained linear system is solved during the inner iteration; when the fracture state no longer changes and the relative change of the normal displacement vector is less than a preset convergence tolerance, it is determined that the double iteration process converges, and the fracture width is output; and the fracture conductivity is evaluated according to the fracture width, so that the key technical bottleneck in large-scale fracture closure calculation can be solved, and a basis for fracturing effect evaluation is provided.

[0016] The FMM acceleration technology adopted in the application reduces the calculation complexity from O(N 2 ) to O(N log N) and the memory requirement from O(N 2 ) to O(N), solves the calculation bottleneck of the traditional method in processing large-scale fracture grids, supports the calculation of more than ten thousand grid units, and provides technical support for the simulation of complex fracture networks. 2

[0017] Further, the double iteration algorithm designed in the application solves the convergence problem in fracture closure calculation, and the strategy of separating the outer state recognition from the inner equation solving guarantees the physical rationality and numerical stability, and is suitable for complex working conditions of fracture state changes.

[0018] Further, the multi-radius spherical proxy function and the grid alignment technology proposed in the application avoid the singularity problem of the traditional plane proxy function, improve the numerical accuracy of the far-field approximation, and reduce the calculation amount by using symmetry.

[0019] Further, compared with the traditional fracture closure analysis method, the application can simulate the fracture closure behavior under complex geometric shapes based on the complete displacement discontinuity method theory, is suitable for the fracturing effect evaluation of various unconventional reservoirs, and has stronger applicability.

[0020] ​Further, the application introduces the projection technology and PCG iteration in the constraint linear system solution, effectively processes the displacement constraint condition of the closed unit, ensures the consistency of the numerical solution and the physical constraint, and is the basis for evaluating the fracture conductivity. BRIEF DESCRIPTION OF DRAWINGS

[0021] Figure 1 A schematic diagram of the evaluation method of hydraulic fracturing effect based on the FMM accelerated iterative algorithm provided by the application is shown in the figure. Figure 2 A comparison chart of the calculation time of the direct method and the FMM method is shown in the figure. Figure 3 A comparison chart of the memory requirement of the direct method and the FMM method is shown in the figure. Figure 4 A chart of the relationship between the acceleration ratio and the number of units of the direct method and the FMM method is shown in the figure. Figure 5 A comparison chart of the fracture width of the direct method and the FMM method is shown in the figure.

[0022] The implementation of the object, functional characteristics and advantages of the application will be further described with reference to the embodiments and the accompanying drawings. DETAILED DESCRIPTION

[0023] In the embodiments of the application, the term "and / or" describes the association relationship of the associated objects, and indicates that there can be three relationships, for example, A and / or B can represent the three cases of A alone, A and B together, and B alone. The character " / " generally represents an "or" relationship between the associated objects before and after it.

[0024] It should be noted that the terms "first", "second" and the like in the specification and claims of the application and the above-mentioned drawings are used to distinguish similar objects, and do not necessarily indicate a specific order or sequence.

[0025] In the embodiments of the application, the term "a plurality of" means two or more, and other quantifiers are similar.

[0026] In order to make the object, technical scheme and advantages of the embodiments of the application more clear, the embodiments of the application will be described in detail below with reference to the drawings. However, those skilled in the art can understand that in the embodiments of the application, many technical details are proposed in order to make the reader better understand the application. However, even without these technical details and various changes and modifications based on the following embodiments, the technical scheme claimed in the application can be implemented. The division of the following embodiments is for the convenience of description, and should not constitute any limitation on the specific implementation mode of the application, and the embodiments can be combined and referred to each other without contradiction.

[0027] The application provides a hydraulic fracturing effect evaluation method based on an FMM accelerated iteration algorithm. Figure 1 As shown in the figure, Figure 1 The flowchart of the hydraulic fracturing effect evaluation method based on the FMM accelerated iteration algorithm is shown.

[0028] At step S100, a three-dimensional fracture grid model is constructed according to fracture geometric parameters and rock mechanical parameters. The fracture geometric parameters include fracture half-length, fracture half-height, spatial step length, etc. The rock mechanical parameters include rock elastic modulus, Poisson's ratio, etc. According to the fracture geometric parameters and the rock mechanical parameters, the fracture surface is discretized in a three-dimensional space to generate a planar grid model containing N units.

[0029] It should be noted that the spatial step length controls the calculation accuracy. The smaller the grid is, the denser the grid is, and the higher the calculation accuracy is. However, the calculation amount also increases.

[0030] At step S200, a calculation equation for calculating the fracture width is constructed based on the displacement discontinuity method, wherein the calculation equation includes a coefficient matrix representing the elastic interaction between the units of the three-dimensional fracture grid model; and the hierarchical compression of the coefficient matrix is performed by using the kernel-independent fast multilevel cell method.

[0031] The control equation of the displacement discontinuity method (DDM) is: .

[0032] A is a coefficient matrix, ; the elements A ij of the coefficient matrix are calculated by a kernel function K ij , which represents the elastic interaction between the units. K ij defines the normal stress at the center point of the source unit i caused by the unit normal displacement discontinuity on the source unit j in an infinite elastic medium, which is the theoretical basis of the displacement discontinuity method.

[0033] D n ∈R N is a normal displacement vector to be solved, i.e., a fracture width; p∈R N is a fluid pressure vector in the fracture; ∈R N is an initial stress vector.

[0034] 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 K ij is: ; wherein Cr=G / (4π(1-ν)); G is a shear modulus. ν is Poisson's ratio; b is the half size of the unit, a = b = Δx / 2; (△x,△y,△z) is the relative position vector between the centers of two units; ; ; J 11 and J 21 denote the basic integral kernel function defined under the mechanical model.

[0035] Generally, there is a calculation bottleneck in directly processing the dense matrix A, and the present application uses the kernel-independent fast multi-level son method for acceleration.

[0036] Specifically, for the near-field interaction in the three-dimensional fracture grid model, the kernel function K ij as the element A ij of the coefficient matrix, the corresponding calculation equation is: ; For the far-field interaction in the three-dimensional fracture grid model, the kernel-independent fast multi-level son method is used to provide the multiplication operator of the implicit matrix vector, and the multiplication operator is used as the calculation equation, and the multiplication operator is ; wherein A is the coefficient matrix, ; D n is the normal displacement vector to be solved, i.e. the fracture width; p is the fluid pressure vector in the fracture; is the initial stress vector; F() is the multiplication operator of the implicit matrix vector; i is the source unit of the three-dimensional fracture grid model, and j is the field unit of the three-dimensional fracture grid model.

[0037] More specifically, an octree spatial partition structure covering all fracture units is constructed, and the number of units contained in the leaf node is controlled by the parameter occ, wherein occ is set to 256, and due to the row and column point merging sorting, the actual number of units contained in each leaf node is about 128; for the near-field interaction (direct neighbor nodes in the tree structure), the matrix elements of the coefficient matrix A are obtained by directly calculating the kernel function; for the far-field interaction, interpolation decomposition is used for low-rank compression: the interaction matrix is calculated by using the proxy point set, the interpolation decomposition is performed to obtain the skeleton points and the redundant points, and a hierarchical compression structure is constructed; through the upward scanning (compression), interaction calculation and downward scanning (refinement) processes, the output of the operator F() is obtained.

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

[0039] The construction method of the proxy point set of the core-independent fast multi-level sub-method includes steps S210 to S240.

[0040] The basic proxy point set P0 is constructed by spherical triangle subdivision at step S210, and a subset of 12 points is randomly sampled from P0 as a basic direction vector. , wherein is a point on the unit sphere. A subset of 12 points is randomly sampled from as a basic direction vector. .

[0041] The radius parameter geometry R j ∈[1.5,2.5], j=1,2, …, k is defined at step S220, and a multi-radius spherical shell proxy point set P multi is constructed by scaling each basic direction vector by the radius value.

[0042] A set of radius parameters {R1, R2, …, R k} is defined, where R j ∈[1.5,2.5], j=1,2, …, k. 12 radius values are selected at equal intervals, and a multi-radius spherical shell proxy point set P multi is constructed by scaling the basic proxy point set , i.e. each basic direction vector by the radius value, to form 144 proxy points (12 basic directions x 12 radii).

[0043] .

[0044] The multi-radius spherical shell proxy point set P multi is aligned with the three-dimensional grid at step S230, and all transformed points form the aligned point set P aligned , ensuring the geometric consistency of the proxy points with the discrete grid. Let △h be the spatial discrete step of the fracture plane. For each point q= (q x , q y , q z ) ∈ P multi in the multi-radius spherical shell proxy point set P multi , a grid alignment transformation is performed, and the transformation is as follows: ; ; .

[0045] Wherein, round() is a rounding function. All transformed points constitute an aligned point set P aligned .

[0046] The proxy point technology of the general FMM is designed for smooth field problems such as electrostatic field. The DDM kernel function involved in hydraulic fracturing has nonlinearity, high gradient and singularity in the near field, and the general method cannot accurately approximate.

[0047] The present application adopts the targeted multi-radius spherical surface sampling: this technology improves the approximation accuracy of the FMM for the DDM kernel function in the near-far field transition zone. The DDM kernel function exhibits nonlinearity and high gradient variation in the near field, and the traditional single spherical proxy point is difficult to accurately capture. Multi-radius spherical shell sampling covers a range of space, avoids singularity, and provides a high-precision interpolation basis. Combining "multi-radius spherical surface sampling" with "grid alignment technology" to optimize the performance of FMM is a solution to the unique problem of near-field accuracy and stability of DDM kernel function.

[0048] The aligned point set P aligned is obtained at step S240 final .

[0049] All duplicate coordinate points are removed from the aligned point set P aligned to form the final FMM proxy point set P final for subsequent fast multipole calculation.

[0050] Multi-radius spherical surface sampling is introduced into the DDM simulation of hydraulic fracturing cracks. By constructing a spherical shell instead of a spherical surface, it systematically improves the approximation accuracy of the kernel function in the entire far field, which is the basis for ensuring the accuracy of crack closure judgment. The grid alignment technology forces the mathematical proxy point set and the discrete grid in the physical space to be geometrically aligned, enhancing the numerical stability and computational efficiency of the algorithm, and is the key to achieving efficient O(Nlog N) calculation.

[0051] At step S300, the rock mechanics parameters are obtained, and the mechanics parameters and the crack state of all elements in the three-dimensional fracture grid model are initialized. The rock mechanics parameters include rock elastic parameters, including elastic modulus E, Poisson's ratio v, and shear modulus; The mechanics parameters include fluid pressure distribution, and the initialization formula is as follows: ; Wherein, p center is the fluid pressure of the central region; i is the source element of the three-dimensional fracture grid model; p edge is the fluid pressure of the edge / boundary part; is a random fluctuation term, ; L center is the lateral extent of the fracture core region; The fracture state state of all elements in the three-dimensional fracture mesh model is a binary state vector, state∈{0,1} N , and the fracture state state of all elements is set to open at the initial time, state(i)=1.

[0052] The fracture width is solved by double iteration of the calculation equation at step S400, wherein the double iteration includes outer iteration and inner iteration. The fracture state is updated during the outer iteration, and the constrained linear system is solved during the inner iteration.

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

[0054] The iteration is performed on each element of the three-dimensional fracture mesh model at step S410, and if the source element i∈I open satisfies , the fracture state state of the source element i is changed from open to closed.

[0055] Let I open and I closed represent the index set of the elements in open and closed states in the current iteration step, respectively. In this embodiment, I m.

[0056] If the source element i∈I closed satisfies , the fracture state state of the source element i is changed from closed to open at step S420; wherein D n (i) is the normal displacement of the source element i; ε d is a preset displacement threshold; is the contact stress of the source element i; p(i) is the fluid pressure; is the initial stress; ε s is a preset stress threshold. In this embodiment, .

[0057] D n (i)=0 when the source element i has closed elements at step S430, and the displacement constraint of the closed elements is handled by a constraint projection function, and the preconditioned conjugate gradient method is used to solve the modified linear system. The definition of the constraint projection function is: Diagonal projection matrix , whose diagonal elements S ii , S ij are determined by the state vector The constraint projection function is taken as an operator : R N → R N , which acts on any input vector to perform the following operation: Where F() is an implicit matrix-vector multiplication operator; state(i)=1 represents that the unit i is in the open state, and state(i)=0 represents that the unit i is in the closed state.

[0058] The operator ensures that the displacement of the closed unit is forced to be zero; the equation corresponding to the closed unit is removed; and the far-field action is accelerated by FMM.

[0059] The projection matrix S is fused with the FMM implicit operator F, so that the displacement constraint processing is realized without explicitly constructing, modifying or storing the dense coefficient matrix A.

[0060] At step S440, the convergence condition is met. The convergence condition is met at the same time as follows: The crack state state changes converge, and ||state （k） -state （k-1） ||=0; the norm of the state vector (the sum of the absolute values of each component) is zero, indicating that the open / closed state of all units no longer changes.

[0061] The normal displacement of the field unit converges relatively, and the formula is as follows: Where k is the iteration number; state （k） is the crack state obtained in the kth iteration; D n (k) is the normal displacement obtained in the kth iteration; tol is a preset convergence tolerance. The relative change of the normal displacement vector is less than the preset convergence tolerance tol, and the maximum iteration number is set to k max =1000, wherein tol is 10 −8 .

[0062] The present application can ensure the most real solution of the crack width through double iteration.​​​​

[0063] At 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, it is determined that the double iteration process converges, and the crack width is output.

[0064] At step S600, the crack conductivity is evaluated according to the crack width. The evaluation of the crack conductivity according to the crack width is a conventional way, which is not specifically introduced here. The larger the crack width, the greater the crack conductivity.

[0065] Example: A three-dimensional crack model is used for numerical verification, and the case parameters are set as follows: the crack half-length is 150 m, the half-height is 30 m, the rock elastic modulus is 30.0 GPa, the Poisson's ratio is 0.2, and the initial ground stress is 10 MPa. The fluid pressure distribution in the crack is non-uniform, the fluid pressure in the central region (|x|≤100 m) is 10 MPa, the fluid pressure in the edge region (|x|>100 m) is 9.5 MPa, and a random disturbance is superimposed to simulate the actual working condition. The spatial step dx is set to 3 m, 2 m, 1 m and 0.5 m, respectively, and the corresponding grid element numbers are 2121, 4681, 18361 and 72721, respectively. The maximum number of iterations is set to 1000, and the convergence tolerance is 10 −8 .

[0066] Solution: According to the evaluation method of hydraulic fracturing effect provided by the FMM accelerated iteration algorithm, first, a crack geometric model is constructed and the mechanical parameters are initialized, then the FMM accelerated displacement discontinuity method is used to construct the coefficient matrix, wherein the proxy function adopts a multi-radius spherical sampling technology, the radius range is [1.5, 2.5], and the number of proxy points p=144. The crack closure problem is solved by a double iteration algorithm, the outer iteration updates the crack state, and the inner iteration solves the constrained linear system by PCG.

[0067] The calculation time comparison is shown in Figure 2 As the number of grid elements increases from 2121 to 72721, the calculation time of the FMM method increases from 0.1 s to 4.1 s, while the calculation time of the direct method is already very long when the number of elements exceeds 4681. At the scale of 18361 elements, the speedup ratio (SR) of the FMM method compared with the direct method reaches 45 times, which reflects the advantage of the present application in large-scale calculation. Figure 4

[0068] The memory requirement comparison is shown in Figure 3 ​As shown, the memory consumption of the FMM method increases approximately linearly with the grid size, from 6.4 MB at 2121 cells to 367.4 MB at 72721 cells; while the memory consumption of the direct method increases quadratically, reaching 175.3 MB at 4681 cells, and the theoretical memory requirement exceeds 2.5 GB at 18361 cells, which exceeds the bearing capacity of a conventional computer.

[0069] The crack width is compared with Figure 5 As shown, to verify the calculation accuracy of the FMM acceleration method, the present application compares the crack width distribution obtained by the direct DDM method and the FMM acceleration method. By extracting the crack width data along the center line of the crack (y=0) for comparative analysis, it is found that the crack width distribution curves calculated by the two methods are basically coincident, which fully proves that the FMM acceleration method improves the calculation efficiency while maintaining the calculation accuracy.

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

[0071] Obviously, the above-described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, those skilled in the art can make other different forms of changes or modifications without making creative efforts, which should all belong to the protection scope of the present application.

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 comprises a coefficient matrix representing the elastic interaction between the units of the three-dimensional fracture grid model; and the coefficient matrix is compressed hierarchically by using a kernel-independent fast multi-level cell method; Obtain rock mechanical parameters, initialize the mechanical parameters and the crack state of all units in the three-dimensional fracture grid model; Solve the calculation equation by double iteration to obtain the fracture width, wherein the double iteration comprises 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 a preset convergence tolerance, it is determined that the double iteration process converges, and the fracture width is output; According to the fracture width, the fracture conductivity is evaluated.

2. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 1, wherein, The hierarchical compression of the coefficient matrix by using the kernel-independent fast multi-level cell method comprises: For the near-field interaction in the three-dimensional fracture mesh model, the kernel function K ij as the element A ij of the coefficient matrix ; 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, ; D n The normal displacement vector to be solved, i.e. the crack width; p is a fluid pressure vector in the fracture; 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.

3. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 2, characterized in that, Considering the interaction between source element i and field element j in the three-dimensional fracture mesh model, the calculation formula of kernel function K ij is as follows: ； Wherein, Cr=G / (4π(1-ν)); G is the shear modulus; ν is the Poisson's ratio; b is the half size of the unit, a=b=△x / 2; (△x,△y,△z) is the relative position vector between the centers of two units; ； ； J 11 With J 21 denotes the basic integral kernel function defined under this mechanical model.

4. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 2, wherein, For the far-field interaction of the three-dimensional fracture grid model, the kernel-independent fast multi-level cell method provides an implicit matrix-vector multiplication operator, and the construction method of the proxy point set of the kernel-independent fast multi-level cell method in the step of using the multiplication operator as the calculation equation comprises: A basic proxy point set P0 is constructed by spherical triangle subdivision, and a subset containing 12 points is randomly sampled from P0 as a basic direction vector; define radius parameter geometry, R j P = {p1, p2, …, pk} is constructed by scaling each basis direction vector by each radius value, j = 1, 2, …, k multi ; For most of the radius of the spherical shell agent point set P multi Using three-dimensional grid alignment transformation, all the transformed points constitute the aligned point set P aligned , ensure the geometric consistency of the proxy point and the discrete grid; Aligned point set P aligned Final proxy point set P formed after deduplication final .

5. 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 the mechanical parameters and the crack state of all units in the three-dimensional fracture grid model, the rock mechanical parameters comprise rock elastic parameters, including the elastic modulus E, the Poisson's ratio v, and the shear modulus; The mechanical parameters comprise a fluid pressure distribution, and the initialization formula is as follows: ； where p center is the fluid pressure in the central region; i is a source unit of the three-dimensional fracture grid model; p edge Fluid pressure at the edge / boundary site; for the random fluctuation term, ; L center limiting the lateral extent of the central region of the fracture; The fracture state state of all elements in the three-dimensional fracture grid model is a binary state vector, state ∈ {0, 1} N , and the fracture state state of all elements is set to open at the initial time, state(i) = 1.

6. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 5, wherein, The step of solving the calculation equation by double iteration to obtain the fracture width comprises: Iterate over each element of the three-dimensional fracture mesh model, if the source element i∈I open satisfies then change the fracture state state of the source element i from open to closed; If the source unit i e I closed satisfies then the crack state state of the source unit i is changed from closed to open. where D n (i) is the normal displacement of source element i; ε d is a predetermined displacement threshold value; Contact stress for source unit i; p(i) is the fluid pressure; initially stressed; ε s is a predetermined stress threshold; 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. Until the convergence condition is met.

7. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 6, characterized in that, In the step of until the convergence condition is met, the convergence condition satisfies the following two conditions at the same time: state changes converge, ||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） state （k） state （k） state （k） state （k） state （k） state D n (k) Normal displacement obtained for the kth iteration; tol is the preset convergence tolerance.

8. The method for evaluating the hydraulic fracturing effect based on the FMM accelerated iterative algorithm according to claim 6, wherein, The definition of the constraint projection function is as follows: Diagonal projection matrix whose diagonal elements S ii , S ij are determined from the state vector ​ The constraint projection function as an operator : R N → R N which acts on an arbitrary input vector by performing 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