A method for evaluating and optimizing fracturing of deep shale fractured reservoirs based on multi-dimensional tortuosity

Abstract

This invention discloses a method for evaluating and optimizing fracturing in deep shale fractured reservoirs based on multidimensional tortuosity, belonging to the field of unconventional oil and gas development and hydraulic fracturing engineering. Based on reservoir geological and construction parameters, this method constructs a fluid-structure interaction numerical model to simulate fracture propagation and stress evolution during fracturing. By identifying the failure modes of tensile and shear fractures, it extracts fracture path length and dynamic aperture, establishing an effective tortuosity evaluation model incorporating fluid dynamics weights. Furthermore, it introduces factors such as fracture zone approach angle and rock cohesion to construct a comprehensive fracture network complexity index including a geological penalty term. With the goal of increasing reservoir stimulation volume and reducing fracture zone blocking effect, it optimizes the fracturing flow rate and fracturing fluid viscosity. This method overcomes the shortcomings of traditional geometric indices in reflecting conductivity, achieving quantitative evaluation of fracture network complexity and closed-loop optimization of fracturing parameters, and has strong engineering applicability.

Description

Technical Field

[0001] This invention relates to the field of unconventional oil and gas field development and hydraulic fracturing engineering technology, specifically to a method for evaluating and optimizing fracturing deep shale fractured reservoirs based on multidimensional tortuosity, belonging to the field of oil and gas field development technology. Background Technology

[0002] Due to their ultra-low permeability and strong heterogeneity, shale gas reservoirs must rely on horizontal well volumetric fracturing technology for commercial development. As shale gas exploration and development in my country moves towards deeper and ultra-deeper reservoirs, the geological structure becomes increasingly complex, and "natural fracture zones" with unidirectional, high local density characteristics have become a significant geological marker of deep reservoirs. Due to the combined influence of factors such as high ground stress difference and fracture zone approach angle, multiple clusters of hydraulic fractures are prone to complex nonlinear mechanical behaviors such as "capture, crossing, turning and blocking" when encountering fracture zones. [Yi Liangping, Yang Changxin, Yang Zhaozhong et al. Influence of natural fracture zones on the propagation of deep shale hydraulic fracturing fractures[J]. Natural Gas Industry, 2022, 42(10): 84-97.] [Zeng Bo, Wang Yufeng, Song Yi et al. Influence of natural fractures in deep shale reservoirs in southern Sichuan Basin on volumetric fracturing fracture network[J]. Daqing Petroleum Geology and Development, 2025, 44(3): 168-174.] [HU YQ, WANG YF, WANG Q, et al. The propagation laws of hydraulic fractures under the influence of natural fracture zones[J]. Physics of Fluids, 2025, 37(1): 016616.]. This intense non-uniform expansion not only makes it difficult to accurately predict reservoir stimulation volume (SRV), but also causes numerous engineering problems such as local fracture network failure and ineffective proppant placement. In order to clarify the evolution law of multi-cluster fractures under the influence of fracture zones, it is first necessary to accurately quantify the complexity of the fracture network. Therefore, constructing a quantitative evaluation index for complex fracture morphology that conforms to the physical meaning of fluid-structure interaction has become an urgent research need for deep shale gas development.

[0003] In current engineering practice and existing research, the reservoir stimulation area and absolute value of total fracture length obtained from microseismic monitoring inversion, or the simple macroscopic geometric tortuosity (i.e., the ratio of total path length to straight-line distance) are often used to evaluate the complexity of fractures and the fracturing effect [Zhang B., Guo T., Chen M., et al. Research on fracture propagation of hydraulic fracturing in a fractured shale reservoir using a novel CDEM-based coupled HM model[J]. Comput. Geotech. 2024, 168, 106-170.]. However, these methods have extremely obvious limitations in practical applications: First, the "geometric length-only" approach is prone to misjudgment. Traditional evaluations only consider the superposition of absolute fracture lengths, without considering the essential differences in rock fracture modes. In actual fracturing, there is an order-of-magnitude difference in the opening and conductivity between the shear fractures activated within the fracture zone and the tensile fractures splitting within the matrix. Superimposing the two with equal weights can easily result in an "artificially high" stimulation volume, leading to serious underestimation of benefits. Secondly, there is a lack of clear hydrodynamic transmission significance. The existing macroscopic tortuosity is essentially only equivalent to the scaling factor of a two-dimensional geometric polygonal line, which seriously deviates from the cubic law of non-Newtonian fluids within the fracture and cannot achieve objective lateral comparisons under different displacements, different fracturing fluid viscosities, and different proppant types. Thirdly, existing conventional numerical simulations (such as finite element model (FEM) and cohesive model (CZM)) cannot scientifically close the loop between the "dynamic evolution evaluation" and "on-site parameter control" of the fracture network, resulting in the inability to achieve real-time process intervention in response to the strong blocking effect of natural fracture zones.

[0004] Furthermore, some cutting-edge studies have proposed using extended finite element method (XFEM) or discrete fracture network (DFN) models to evaluate the fracturing degree of fracture zones. However, such evaluation systems suffer from three significant technical bottlenecks in practical industrial applications: First, their numerical algorithms are prone to mesh distortion and computational non-convergence when dealing with numerous bifurcations and intersections of multi-cluster fractures with tens of millions of meshes, resulting in severely distorted evaluation results and low efficiency. Second, their evaluation models employ static geological assumptions and lack dynamic penalty mechanisms for extreme values ​​of natural fracture approach angles (such as the 45° optimal window) and rock cohesion, limiting the method to post-event review evaluations and failing to realize the engineering value of real-time online control of fracturing parameters (displacement, viscosity). Third, a single fracture length evaluation cannot distinguish between "effective fluid sweep and support" and "pure shear slip without conductivity" from a physical and mechanical perspective, easily leading to overestimation of the modification effect in complex deep geostress environments.

[0005] In summary, existing quantitative evaluation methods for the complexity of hydraulic fractures in natural fracture zones all have certain limitations, and a dynamic evaluation system capable of accurately identifying fracture modes and constructing effective hydrodynamic transmission significance is still lacking. Therefore, in the context of large-scale application of multi-cluster fracturing in deep shale gas, it is imperative to construct a reasonable index for quantifying the complexity of fracture networks. To accurately characterize the effective conductivity of fracture networks and effectively overcome the blocking effect of natural fracture zones, it is urgent to establish a fracturing evaluation and intelligent optimization method with clear physical meaning and integrating multi-dimensional dimensionality reduction algorithms to fill the gap in oilfield development technology in this field. Summary of the Invention

[0006] This invention primarily overcomes the shortcomings of existing technologies by proposing a multidimensional tortuosity-based evaluation and optimization method for hydraulic fracturing in deep shale fractured reservoirs. This method considers the combined effects of fluid flow, formation stress changes, local fracture aperture, fracture mode differences, and the approach angle of natural fracture zones during fracturing. Based on elasticity, seepage mechanics, numerical simulation principles, and fluid-structure interaction theory, a model for evaluating the complexity of hydraulic fracturing in deep shale fracture zones is established. Using this model, the complexity of the comprehensive fracture network can be predicted based on construction and geological parameters, and a scientifically reasonable fracturing displacement and fracturing fluid viscosity can be determined. This proposed method fills a gap in existing research theories and technologies.

[0007] The technical solution provided by this invention to solve the above-mentioned technical problems is: a method for evaluating and optimizing the fracturing of deep shale fractured reservoirs based on multidimensional tortuosity, comprising the following steps ( Figure 1 ): S1. Based on the geological research of deep shale gas, obtain the geological characteristic parameters of the target shale gas, including formation pressure, matrix porosity, permeability, natural fracture zone parameters (approach angle, bandwidth, distribution density of natural fractures) and rock mass mechanical parameters, and determine the initial conditions and boundary conditions. S2. Based on the above geological characteristic parameters of deep shale gas, a seepage-geomechanical coupled model is established by combining the continuity equation of the finite-discrete element method (FDEM) with the fluid dynamics control equation, and the distribution of deformation, fracturing and pore pressure in the formation during the fracturing process is obtained. S3. Based on the displacement and stress field distribution obtained by the coupling model, the failure modes of joint units in deep shale gas reservoirs after fracturing are extracted, and the fracture trajectories, cumulative lengths and real-time fracture widths of tensile failure (matrix cracking) and shear failure (activation of natural fractures) are separated and obtained. S4. Based on the flow-guiding force of fluid within multiple clusters of cracks, an effective tortuosity evaluation model incorporating fluid dynamics weighting factors is derived and established. S5. Based on the geometric characteristics and stress disturbance of natural crack zones, establish a comprehensive crack network complexity index prediction model that includes a crack zone induced deviation penalty term. S6. Based on the comprehensive fracture network complexity index calculated in steps S4 and S5, with the optimization objective of maximizing the modified volume and reducing the blocking effect of the fracture zone, the optimal combination of single cluster displacement and fracturing fluid viscosity during the fracturing process is obtained.

[0008] A further technical solution is that the specific process of step S1 is as follows: based on the geological research of deep shale gas, obtain the existing geological characteristic parameters of the target shale gas, and determine the initial conditions and boundary conditions. The specific parameters include formation pressure, matrix and natural fracture cohesion, tensile strength, construction discharge rate, total liquid volume, liquid viscosity, number of fracture clusters, etc.

[0009] A further technical solution is that the specific process of step S2 is as follows: S201, Discretize the target reservoir into triangular units and zero-thickness four-node joint units ( Figure 2 Combining Newton's second law with the potential function and contact force, a global governing equation is established to obtain the formation deformation distribution: ; in, m For node quality, N·s 2 / m; u Let the displacement be m; C The damping coefficient is N·s / m; f c The contact force between discrete elements is expressed in Pa. f d The deformation force of the triangular element is Pa; f j The joint unit bonding force is expressed in Pa. f ext Let Pa be the external force vector of fluid pressure. S202. The relationship between flow rate and crack opening is updated in real time using the cubic law: ; in, q For fluid flow rate, m 2 / s; a Let m be the average aperture of the joint element; The viscosity of the liquid is Pa·s; The pressure difference at the nodes is expressed in Pa. L f The node spacing is in meters (m). A further technical solution is that the specific process of step S3 is as follows: the specific judgment logic for tensile failure and shear failure is as follows: when the normal tensile stress on the joint unit reaches the local tensile strength of the rock and the normal opening displacement reaches the critical value, the system determines and records it as tensile failure; when the tangential stress on the joint unit reaches the local shear strength of the rock and the tangential slip displacement reaches the critical value, the system determines and records it as shear failure; if both occur simultaneously, it is recorded as a tensile-shear composite failure, and it is split and counted according to the displacement ratio. L s and L t middle( Figure 3 ); A further technical solution is that the specific process of step S4 is as follows: breaking through the traditional method of relying solely on geometric length calculation, comprehensively extracting the cumulative length of shear cracks. L s Cumulative length of tensile cracks L t And its average opening, deriving the calculation of effective weighted tortuosity. T e ( Figure 4 ): ; in, T e For effective weighted tortuosity, dimensionless; L s The cumulative length of the shear failure crack, in meters; L t The cumulative length of the tensile failure crack, in meters; l The maximum projected distance, in meters, spanned by multiple clusters of hydraulic fractures in the direction perpendicular to the horizontal wellbore. and Here, represents the average aperture of shear and tensile cracks, respectively, in meters (m). , The effective conductivity weighting coefficient, determined by the type of fracturing fluid and proppant, is dimensionless. a 0 is the reference baseline opening, which is usually taken as the initial average closure opening of the natural fractures in the target reservoir, or directly as a unit constant with uniform dimensions. The purpose is to convert the dynamic evolution of the real fracture opening into a dimensionless multiple relative to the baseline state, thereby ensuring that the comprehensive tortuosity has the ability to make horizontal comparisons and evaluations across scales and blocks.

[0010] A further technical solution is that the specific process of step S5 is as follows: S501. Define a fracture zone induced deviation penalty term to characterize the effective energy blocking caused by large approach angles and high cohesion: ; in, The actual approach angle of the fracture zone is the angle between the natural fracture and the vertical wellbore, expressed in °. The optimal approximation angle constant (set to 45°) is °; C nf and C matrix These represent the cohesion between the natural fracture zone and the matrix, in MPa. S502. Combining the inter-cluster interference theory of multi-cluster fracturing, a logarithmic function of topological intersection gain is defined; coupling these two functions with the horizontal principal stress difference yields the final comprehensive fracture network complexity exponential evolution equation: ; in, C fni The overall complexity index of the seam mesh is dimensionless. N cross The total number of intersections of multiple crack clusters; N main The number of main pier clusters; cluster; The difference in horizontal principal stresses is expressed in MPa. The maximum horizontal principal stress is expressed in MPa. This is a penalty term for deviation induced by cracks; This formula achieves for the first time an explicit nonlinear characterization of the fracture zone approximation characteristics (with 45° as the extreme singularity) and the macroscopic modification volume efficiency; it also integrates the fracture network complexity index. C fni The curve showing the change in the approach angle of the natural fracture zone is as follows: Figure 5 As shown, this diagram visually reveals C fni The nonlinear decay law that occurs as the approximation angle deviates from 45°.

[0011] A further technical solution is that the specific process of step S6 is as follows: when the model predicts... C fni When the pressure falls below the set engineering economic threshold (indicating high cohesion or large approach angle obstruction by the crack zone), the optimization module automatically triggers the compensation mechanism, outputting a closed-loop feedback command to increase the single-cluster discharge rate or increase the liquid viscosity to improve the net fluid pressure. Figure 6 , Figure 7 ).

[0012] The beneficial effects of this invention are: using this method to predict the comprehensive tortuosity evolution of fracture networks in deep shale gas reservoirs during multi-cluster fracturing, and establishing a fusion... T e and penalty items The exponential model accurately identifies advantageous modification windows with an approximation angle close to 45° and low cohesion. This proposed method fills the gap in existing theoretical and technical research on post-fracturing control of deep fracture zones and provides a reliable analytical research method for fracturing optimization of complex reservoirs. Attached Figure Description

[0013] Figure 1 Flowchart of evaluation and optimization method for fracturing deep shale fractured reservoirs based on multidimensional tortuosity; Figure 2 Construction diagram of finite discrete element numerical computation unit; Figure 3 Comparison of joint element failure modes at different approximation angles; Figure 4 Schematic diagram for calculating the evaluation index (torsion) of multiple clusters of hydraulic fractures; Figure 5 Comprehensive stitching complexity index C fni Curve showing the change in approximation angle; Figure 6 Response surface diagram of the correlation between optimal construction discharge rate and geological parameters; Figure 7 Response surface plot of optimal fracturing fluid viscosity in relation to geological parameters. Detailed Implementation

[0014] To make the objectives, technical solutions, and advantages of the present invention clearer, the present invention will be further illustrated below with reference to the accompanying drawings and embodiments.

[0015] This invention provides a method for evaluating and optimizing fracturing in deep shale fractured reservoirs based on multidimensional tortuosity. Its core lies in using a finite-discrete unit bottom-layer physical framework and fluid dynamics-weighted dimensionality reduction to handle the non-uniform evolution of multiple fracture clusters within the fracture zone, and efficiently embedding it into a dynamic evaluation framework with real-time feedback of fracturing parameters. The specific implementation steps are as follows (…). Figure 1 ): S1. Based on the geological and engineering characteristics of the target reservoir, obtain basic mechanical data of deep shale and detailed parameters for natural fracture zones; S2. Based on the finite-discrete element method, a seepage-geomechanical coupled model is established to calculate the deformation, fracturing and pore pressure distribution in the formation during the fracturing process; S3. Based on the fluid-structure interaction seepage mechanics model, extract the tensile and shear failure trajectories of multiple crack clusters and the real-time dynamic mesh opening; S4. Establish an effective weighted dynamic evaluation model for tortuosity that considers fluid dynamics transmission and geometric topology; S5. Construct a comprehensive fracture network complexity index tensor that couples geological penalties of natural fracture zones with stress interference between clusters; S6. Closed-loop control and dynamic solution of fracturing strategy based on multidimensional threshold matrix; In this invention, the specific acquisition requirements for step S1 are as follows: Obtain the geological characteristic parameters and fracturing operation parameters of the existing target unconventional oil and gas reservoir. Specific parameters include: matrix rock elastic modulus, Poisson's ratio, maximum / minimum horizontal principal stress, natural fracture zone approach angle, matrix-fracture cohesion, fracturing fluid system viscosity, operation displacement, number of fracturing clusters per segment, and response projection reference length. Simultaneously, import bedding and microfracture distribution characteristics using modeling software. Figure 2 ).

[0016] In this invention, the specific process of step S2 is as follows: Traditional crack morphology monitoring usually only extracts the length of a single macroscopic envelope. This invention innovatively introduces the finite discrete element (FDEM) underlying motion differential calculus; tracks the displacement field and stress field of the target mesh node, accurately captures the critical singular points of element opening and slip, and uses the cubic law to update the dynamic response relationship between flow rate and crack opening of non-Newtonian fluid in real time.

[0017] In this invention, the specific process of step S3 is as follows: separating the tensile failure path within the matrix from the shear failure path within the natural fracture zone ( Figure 3 The average aperture of tensile cracks and the average aperture of shear cracks at any given time are dynamically extracted. These two parameters not only reflect the apparent crack network volume, but are also the core boundary conditions for constructing effective flow conduction weights.

[0018] In this invention, the specific process of step S4 is as follows: S401. In order to quantitatively determine the true hydrodynamic physical scale of multiple crack clusters, this invention breaks through the traditional macroscopic ratio at the purely geometric level and substitutes the classification path length and real-time opening extracted in S3 into the fluid energy equivalent framework. S402. For any complex interwoven mesh morphology, its effective transmission and transformation capabilities depend not only on the total path length but also on the strict constraints of the mesh width and failure mode. Based on this, a dynamic calculation model for the effective tortuosity weighted by fluid dynamics is derived. Figure 4 ): ; in, T e For effective weighted tortuosity, dimensionless; L s The cumulative length of the shear failure crack, in meters; L t The cumulative length of the tensile failure crack, in meters; lThe maximum projected distance, in meters, spanned by multiple clusters of hydraulic fractures in the direction perpendicular to the horizontal wellbore. and Here, represents the average aperture of shear and tensile cracks, respectively, in meters (m). , The effective conductivity weighting coefficient, determined by the type of fracturing fluid and proppant, is dimensionless. a 0 is the reference baseline opening, which is usually taken as the initial average closure opening of the natural fractures in the target reservoir, or directly as a unit constant with uniform dimensions. The purpose is to convert the dynamic evolution of the real fracture opening into a dimensionless multiple relative to the baseline state, thereby ensuring that the comprehensive tortuosity has the ability to make horizontal comparisons and evaluations across scales and blocks.

[0019] In this invention, the specific process of step S5 is as follows: S501. Traditional evaluation models typically neglect the nonlinear hindering effect of geostress deflection and fracture zone properties on fracture penetration. This invention, based on weighted effective tortuosity, innovatively introduces a local geological penalty and spatial topological gain mechanism; it defines a fracture zone-induced deviation penalty term to characterize the effective energy blockage caused by large approach angles and high cohesion. ; in, The actual approach angle of the fracture zone is the angle between the natural fracture and the vertical wellbore, expressed in °. The optimal approximation angle constant (set to 45°) is °; C nf and C matrix These represent the cohesion between the natural fracture zone and the matrix, in MPa. S502. Combining the inter-cluster interference theory of multi-cluster fracturing, a logarithmic function of topological intersection gain is defined. Coupled with the difference in horizontal principal stress, the final comprehensive fracture network complexity exponential evolution equation is obtained: ; in, C fni The overall complexity index of the seam mesh is dimensionless. N cross The total number of intersections of multiple crack clusters; N main The number of main pier clusters; cluster; The difference in horizontal principal stresses is expressed in MPa. The maximum horizontal principal stress is expressed in MPa. This is a penalty term for deviation induced by cracks; This index is the first to achieve an explicit nonlinear characterization of the approximation characteristics of natural fracture zones and the macroscopic volumetric efficiency of alteration; such as Figure 5As shown, as the approach angle of the natural fracture zone increases from 20° to 80°, the overall fracture network complexity index increases. C fni It exhibits an inverted V-shaped trend of first increasing and then decreasing, reaching a maximum value near 45°; this curve demonstrates from a physical and mechanical perspective that when the approach angle is extremely small or extremely large, hydraulic fractures are prone to single shear slip or direct penetration, leading to increased fracture network complexity. C fni The 45° approach angle can induce the most complex tension-shear composite fracture, which is the optimal transformation window for deep shale gas fracturing construction.

[0020] In this invention, step S6 specifically involves: constructing a multi-dimensional decision space in the form of displacement within the framework of full life-cycle fracturing simulation and real-time construction control; and employing a sequential iteration and dynamic evaluation mechanism within each fracturing construction time step. S601. First, input real-time geological logging data and pressure data, and calculate the current time. C fni The system automatically checks whether these two indicators meet the maximum convergence tolerance for reservoir stimulation volume, along with the extreme values ​​of the penalty term. S602, if C fni If the temperature continues to rise steadily and deviates very little from the penalty term (fracture zone approach angle approaches 45°), it is identified as a "fracture zone induced enhancement zone," and the system outputs a command to maintain the design flow rate and low-viscosity fracturing fluid to fully activate shear slip; if C fni When the boundary becomes stable, it is determined to be a steady-state matrix expansion region. At this point, stable interference occurs between clusters, macroscopic cracks expand at equal lengths, and the process proceeds smoothly according to the original design parameters. The node bifurcation rate is also monitored more frequently. S603, If detected within the time step... C fni If an exponential drop occurs and the penalty term approaches the peak boundary (large angle or strong cohesion conditions), it is determined that the system has entered the "strong blocking malignant penetration zone". The system triggers a real-time intervention mechanism and issues commands to increase the injection rate of a single cluster and switch to high-viscosity hydroplaning for forced pressure increase until the crack crosses the crack zone again and meets the complexity exponential convergence condition.

[0021] Based on the fluid-structure interaction simulation and multidimensional tortuosity quantitative evaluation steps described in this invention, firstly, by accurately capturing the complex mechanical evolution behavior of multiple clusters of hydraulic fractures in deep natural fracture zones—including their capture, crossing, and obstruction—and the underlying fluid-structure interaction physical mechanism, a solid theoretical foundation is provided for scientifically assessing the fracture zone induction effect and accurately predicting the complexity of the fracture network for effective reservoir stimulation volume (SRV). Secondly, by using fluid dynamics weighting and spatial topology penalty mechanisms (such as approximation angle extreme value constraints) to mathematically reconstruct and reduce the dimensionality of the comprehensive evaluation model of the fracture network, the technical and computational bottleneck of relying solely on the total geometric fracture length to qualitatively judge the stimulation effect at the fracturing site, which easily masks the distribution of ineffective tensile / shear fractures, is completely overcome. Finally, this invention not only significantly improves the computational efficiency and accuracy of field feedback in evaluating fracture morphology evolution under complex geological conditions, but also provides highly valuable technical guidance for the efficient development of deep shale gas fracture zones, the intelligent closed-loop control of fracturing construction parameters (such as displacement and viscosity), and the precise quantification of the final single-well recovery rate (EUR).

[0022] Example 1: Geological data of a deep shale gas well in the southern Sichuan Basin was obtained through field logging. The fracturing process parameters for this well relied heavily on experience and lacked scientifically sound technical guidance. This well is typical of deep shale gas reservoirs, characterized by a burial depth exceeding 3500m, well-developed natural fracture zones, and high formation pressure. This well serves as a good demonstration case for this method. The specific simulation steps are as follows (…). Figure 1 ): 1. Parameter Acquisition: The geological parameters of the reservoir were obtained through field logging: Poisson's ratio of the rock is 0.2, the original formation pressure of the reservoir is 80 MPa, the elastic modulus is 38 GPa, the horizontal principal stress difference is 12 MPa; the fracture zone thickness is 20 m, the natural fracture approach angle is 65°, the matrix cohesion is 20 MPa, and the fracture zone cohesion is 5 MPa; the planned injection rate is 2 m³ / min per cluster, and the fracturing fluid viscosity is 5 mPa·s.

[0023] 2. Construct a geological physical model based on known parameters. Figure 2 Using the finite discrete element fluid-structure interaction differential equation of this invention, discrete difference calculations were performed on the model to obtain the dynamic propagation displacement and real-time opening of the crack during the injection process. a .

[0024] 3. Extract the calculated data and separate the trajectories of tensile failure and shear failure. Figure 3 Substituting into the weighted tortuosity equation in the claims, the calculation is performed. T e ( Figure 4 Since the fracture zone approach angle in this well is 65°, the calculation system identifies a significant fluid blocking effect, which is determined through equations. Activating the penalty factor leads to the final calculatedC fni Below the set effective modification threshold of 2.5 ( Figure 5 ).

[0025] 4. Optimize the feedback system by automatically performing regression calculations based on the objective function. It is recommended to increase the single-cluster fracturing flow rate from 2 m³ / min to 3.5 m³ / min and switch the fracturing fluid system to a linear gel with a viscosity of 30 mPa·s. The combined high flow rate and high viscosity will generate strong net fluid pressure, forcing the two blocked edge fracture clusters to break through the fracture zone barrier and re-energize the simulated output. C fni The area that meets the standards ( Figure 6 , Figure 7 ).

[0026] 5. Following the simulation steps described in this invention, parameters such as formation pressure, velocity distribution, and multidimensional tortuosity evaluation during deep shale gas fracturing can be obtained. These parameters can then be used to intelligently output the optimal combination of fracturing displacement and viscosity, providing closed-loop guidance for on-site construction. Based on the quantitative evaluation and optimization method described in this invention, firstly, by accurately capturing the complex evolutionary behavior of multiple clusters of hydraulic fractures in deep natural fracture zones—"capturing, crossing, and blocking"—and the underlying fluid-structure interaction physical mechanism, a solid theoretical foundation is provided for scientifically assessing the fracture zone induction effect and accurately predicting the complexity of the fracture network in the effective reservoir stimulation volume (SRV). Secondly, fluid is innovatively introduced... The mechanical weighting factor (effective tortuosity) and the fracture zone induced deviation penalty term are used to perform multidimensional nonlinear topological reconstruction on the traditional single static geometric length evaluation model. This completely breaks through the technical bottleneck of relying solely on the total macroscopic fracture length to qualitatively judge the transformation effect at the fracturing site, which easily masks ineffective fractures and leads to false "overestimation" misjudgments. Ultimately, this invention not only significantly improves the accuracy and real-time intervention capability of fracture evolution evaluation under complex geological conditions, but also provides highly valuable technical guidance for the efficient development of deep shale gas fracture zones, the intelligent closed-loop control of fracturing construction parameters (such as displacement and viscosity), and the maximization of the final single-well recovery rate (EUR).

[0027] The above description is not intended to limit the present invention in any way. Although the present invention has been disclosed through the above embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some changes or modifications to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A method for evaluating and optimizing fracturing in deep shale fractured reservoirs based on multidimensional tortuosity, characterized in that, Includes the following steps: S1. Based on the geological research of deep shale gas, obtain the geological characteristic parameters of the target shale gas, including formation pressure, matrix porosity, permeability, natural fracture zone parameters (approach angle, bandwidth, distribution density) and rock mass mechanical parameters, and determine the initial conditions and boundary conditions. S2. Based on the above geological characteristic parameters of deep shale gas, a seepage-geomechanical coupled model is established by combining the continuity equation of the finite-discrete element method (FDEM) with the fluid dynamics control equation, and the distribution of deformation, fracturing and pore pressure in the formation during the fracturing process is obtained. S3. Based on the displacement and stress field distribution obtained by the coupling model, the failure modes of joint units in deep shale gas reservoirs after fracturing are extracted, and the fracture trajectories, cumulative lengths and real-time fracture widths of tensile failure (matrix cracking) and shear failure (activation of natural fractures) are separated and obtained. S4. Based on the flow-guiding force of fluid within multiple clusters of cracks, an effective tortuosity evaluation model incorporating fluid dynamics weighting factors is derived and established. S5. Based on the geometric characteristics and stress disturbance of natural crack zones, establish a comprehensive crack network complexity index prediction model that includes a crack zone induced deviation penalty term. S6. Based on the comprehensive fracture network complexity index calculated in steps S4 and S5, with the optimization objective of maximizing the modified volume and reducing the blocking effect of the fracture zone, the optimal combination of single cluster displacement and fracturing fluid viscosity during the fracturing process is obtained.

2. The method for evaluating and optimizing fracturing deep shale fractured reservoirs based on multidimensional tortuosity according to claim 1, characterized in that, The specific process of step S1 is as follows: Based on the geological research of deep shale gas, obtain the existing geological characteristic parameters of the target shale gas, and determine the initial conditions and boundary conditions. The specific parameters include formation pressure, cohesion between matrix and natural fractures, tensile strength, construction discharge rate, total liquid volume, liquid viscosity, number of fracture clusters, etc.

3. The method according to claim 1, characterized in that, The specific steps of S2 include: S201. Discretize the target reservoir into triangular elements and zero-thickness four-node joint elements. Combine Newton's second law and potential function contact force to establish global control equations and obtain the formation deformation distribution: ； in, m For node quality, N·s 2 / m; Let the displacement be m; C The damping coefficient is N·s / m; f c The contact force between discrete elements is expressed in Pa. f d The deformation force of the triangular element is Pa; f j The joint unit bonding force is expressed in Pa. f ext Let Pa be the external force vector of fluid pressure. S202. The relationship between flow rate and crack opening is updated in real time using the cubic law: ； in, q For fluid flow rate, m 2 / s; Let m be the average aperture of the joint element; The viscosity of the liquid is Pa·s; The pressure difference at the nodes is expressed in Pa. L f The node spacing is in meters (m).

4. The method according to claim 1, characterized in that, The specific steps of S3 include: obtaining the specific judgment logic for tensile failure and shear failure as follows: when the normal tensile stress on the joint unit reaches the local tensile strength of the rock and the normal opening displacement reaches the critical value, the system judges and records it as tensile failure; when the tangential stress on the joint unit reaches the local shear strength of the rock and the tangential slip displacement reaches the critical value, the system judges and records it as shear failure; if both occur simultaneously, it is recorded as a tensile-shear composite failure, and it is split and counted according to the displacement ratio. L s and L t middle.

5. The method according to claim 1, characterized in that, The specific steps of S4 include: comprehensively extracting the cumulative length of shear cracks. L s Cumulative length of tensile cracks L t And its average opening, deriving the calculation of effective weighted tortuosity. T e : ； in, T e For effective weighted tortuosity, dimensionless; L s The cumulative length of the shear failure crack, in meters; L t The cumulative length of the tensile failure crack, in meters; l The maximum projected distance, in meters, spanned by multiple clusters of hydraulic fractures in the direction perpendicular to the horizontal wellbore. and Here, represents the average aperture of shear and tensile cracks, respectively, in meters (m). , The effective conductivity weighting coefficient, determined by the type of fracturing fluid and proppant, is dimensionless. a 0 represents the reference opening, in meters.

6. The method according to claim 1, characterized in that, The specific steps of S5 include: S501. Define a fracture zone induced deviation penalty term to characterize the effective energy blocking caused by large approach angles and high cohesion: ； in, The actual approach angle of the fracture zone is the angle between the natural fracture and the vertical wellbore, expressed in °. The optimal approximation angle constant (set to 45°) is °; C nf and C matrix These represent the cohesion between the natural fracture zone and the matrix, in MPa. S502. Combining the inter-cluster interference theory of multi-cluster fracturing, a logarithmic function of topological intersection gain is defined; coupling these two functions with the horizontal principal stress difference yields the final comprehensive fracture network complexity exponential evolution equation: ； in, C fni The overall complexity index of the seam mesh is dimensionless. N cross The total number of intersections of multiple crack clusters; N main The number of main pier clusters; cluster; The difference in horizontal principal stresses is expressed in MPa. The maximum horizontal principal stress is expressed in MPa. This is a penalty term for deviation induced by crack zones.

7. The method according to claim 1, characterized in that, The specific steps of step S6 include: when the model predicts... C fni When the pressure falls below the set engineering economic threshold, the optimization module automatically triggers the compensation mechanism, outputting a closed-loop feedback command to increase the single-cluster discharge rate or increase the liquid viscosity to improve the fluid net pressure.

Images