Natural fracture development shale reservoir fracturing horizontal shaft casing variable calculation method
By establishing a stress model of the formation-natural fracture-casing system and taking into account the changes in in-situ stress caused by hydraulic fracturing, the DDM method is used to simulate the shear slip length, which solves the problem of the existing technology's failure to accurately predict shale gas horizontal wellbores. To solve the problem of arbitrage risk, more accurate arbitrage prediction and rapid calculation are achieved.
Patent Information
- Application Number
- CN202411504570.4
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2024-10-26
- Publication Date
- 2025-01-24
AI Technical Summary
The prior art fails to consider the ground stress changes caused by hydraulic fracturing when calculating the sheathing of shale gas horizontal wellbore, which leads to the inability to accurately predict the risk of sheathing.
The stress model of the formation-natural fracture-casing system before fracturing construction was established, and the stress field changes caused by formation fracturing were taken into account. The DDM method was used to simulate the shear slip length of the natural fracture zone and calculate the casing variables.
Provide a prediction method that is more in line with the actual construction situation, which can accurately predict the risk of the set change, reduce the dimension and quickly calculate the sliding amount of cutting.
Smart Images

Figure CN119358452A_ABST
Abstract
Description
Technical Field
[0001] The invention belongs to the field of petroleum engineering, and in particular relates to a method for calculating casing variables of a horizontal wellbore in a shale reservoir with natural fractures. Background Art
[0002] Deep shale gas is the main battlefield for natural gas extraction in my country in the future. Horizontal well segmented multi-cluster volume fracturing technology has become the key to its effective development. However, the casing deformation problem (referred to as casing deformation) seriously affects the integrity of the wellbore and the lowering and use of downhole tools, and even leads to the abandonment of reservoir transformation of some fracturing sections, which restricts the large-scale effective development of deep shale gas.
[0003] Existing studies believe that the main factors that induce casing change in shale gas horizontal well transformation are natural fracture slip shearing wellbore, wellbore pipe wall yielding and collapse, and thermal stress damage. Regarding the research on casing change induced by natural fracture slip shearing wellbore, although there is a basic theoretical understanding of the stress failure mechanism of fracture slip and casing change during fracturing and risk prediction, the casing change calculation models established in the current research have not taken into account the impact of the ground stress change caused by hydraulic fracturing, resulting in the inability to accurately predict the risk of casing change when using such models in field engineering design. Summary of the invention
[0004] The object of the present invention is to provide a method for calculating casing variables of a horizontal wellbore in a shale reservoir with natural fractures.
[0005] In order to achieve the above object, the technical solution adopted by the present invention is:
[0006] A method for calculating casing variables of a horizontal wellbore for fracturing a shale reservoir with natural fractures, characterized by comprising the following steps:
[0007] 1) Establish a stress model of the formation-natural fracture-casing system before fracturing construction;
[0008] 2) Establish the stress model for the case where the crack is not fully opened and the stress model for the case where the crack is fully opened;
[0009] 3) Considering the stress field changes caused by formation fracturing, the corresponding pressure increment model is established;
[0010] 4) The DDM method is used to simulate the shear slip length of the natural fracture zone and perform set variable calculations.
[0011] Preferably, the stress model of the formation-natural fracture-casing system in step 1) is:
[0012]
[0013] Where: x , τ xyis the normal stress and shear stress of natural fractures, MPa; σ H , σ h are the maximum and minimum horizontal principal stresses, MPa; θ is the angle between the natural fracture and the maximum horizontal principal stress, rad.
[0014] Preferably, the stress model of the crack not fully opened in step 2) is:
[0015] τ c =(τ xy A f -f max ) / A c (2)
[0016] f max =μσ n A f (3)
[0017] Where: τ c is the shear stress on the casing, MPa; τ xy is the shear stress on the natural fracture, MPa; A f is the natural fracture area, m 2 ;f max is the friction force, N; A c is the cross-sectional area of the intersection of casing and natural fracture, m 2 ; μ is the friction coefficient of the crack surface; σ n is the rock contact normal stress on the fracture surface, MPa.
[0018] Preferably, the stress model for the crack to be fully opened in step 2) is:
[0019]
[0020] The cross-sectional area under the complete cement sheath is:
[0021]
[0022] The cross-sectional area of cement sheath failure is:
[0023]
[0024] Where: τ c is the shear stress on the casing, MPa; τ xy is the shear stress on the natural fracture, MPa; A f is the natural fracture area, m 2 ; A c is the cross-sectional area of the intersection of casing and natural fracture, m 2 ; R is the wellbore radius, m; α is the horizontal wellbore approach angle, °; r 1is the inner diameter of the casing, m; r is the outer diameter of the casing, m; θ is the angle between the natural fracture and the maximum horizontal principal stress, rad.
[0025] Preferably, the step 3) further comprises:
[0026] Establish pore pressure increment model;
[0027]
[0028] The critical pore pressure increment model for natural fracture activation and instability is:
[0029]
[0030] p(x,y)≥p 临界 (9)
[0031] Where: p(x,y) is the gap pressure increment at the plane (x,y); p fi,j (x, y) is the pore pressure generated by the jth segment of the i-th perforation cluster at the plane (x, y); p e is the original formation pressure, MPa; σ H , σ h are the maximum and minimum principal stresses, MPa.
[0032] The stress increment model of the change of three-dimensional ground stress due to the increase of injection pore pressure is:
[0033]
[0034] Δσ xy =σ 诱导mxy (12)
[0035] Where: ν is the Poisson's ratio of the reservoir rock, dimensionless; α is the Biot poroelastic coefficient, dimensionless; p p is the current formation pressure, MPa; σ 诱导mx , σ 诱导my , σ 诱导mxy are the normal stress and shear stress components generated by the crack in the x and y coordinates, MPa.
[0036] Preferably, in step 4), the DDM method is used to simulate the shear slip length of the natural fracture zone, and the corresponding model is:
[0037]
[0038] Where: is the shear stress and normal stress in the i-th fracture unit, MPa; p f is the fluid pressure in the natural fracture zone, MPa; β i is the angle between unit i and the positive x-axis, in rad; is the shear stress and normal stress generated by hydraulic fracture, MPa; is the shear stress caused by fracturing fluid loss, MPa.
[0039] On the other hand, the present invention provides a computing device for calculating horizontal wellbore casing variables for fracturing shale reservoirs with natural fractures, characterized in that it includes: at least one processor and a memory communicatively connected to the processor; the memory stores instructions executable by the processor, and the instructions are executed by at least one processor so that the at least one processor can execute any of the above-mentioned computing methods.
[0040] On the other hand, the present invention provides a computer-readable storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the steps of the computing method as described in any one of the above items.
[0041] In summary, the present invention proposes a method for calculating casing variation of a horizontal wellbore in a shale reservoir with natural fractures. The method of the present invention overcomes the defect that the calculation models in the prior art do not take into account the changes in ground stress caused by hydraulic fracturing, and proposes a casing variation prediction method that is more in line with actual construction conditions, and has the advantages of reducing dimensionality and quickly calculating shear slip. BRIEF DESCRIPTION OF THE DRAWINGS
[0042] Figure 1 It is a schematic diagram of the method flow of the present invention;
[0043] Figure 2 is a graph showing the relationship between the natural fracture area and the wellbore shear stress in an embodiment of the present invention;
[0044] Figure 3 is the relationship diagram between the approach angle of natural fractures and the wellbore shear stress;
[0045] Figure 4 The relationship between the natural fracture approach angle and the wellbore shear stress at different wellbore approach angles;
[0046] Figure 5 The relationship between the approach angle of natural fractures and the wellbore shear stress under different fluid pressures in the fractures;
[0047] Figure 6 It is the relationship diagram between the wellbore position coordinates and the slip amount under different shear stresses;
[0048] Figure 7 The relationship between the wellbore position coordinates and the slip amount under different Young's moduli;
[0049] Figure 8 The relationship between the wellbore position coordinates and the slip amount under different Poisson's ratios;
[0050] Fig. 9 This is an analysis chart of the main controlling factors that affect the set change. DETAILED DESCRIPTION
[0051] The present invention will be described in detail below in conjunction with the accompanying drawings.
[0052] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.
[0053] In the prior art, the change of ground stress caused by hydraulic fracturing is ignored when predicting casing variation, which is inconsistent with the actual physical process. Based on this, the present invention provides a method for calculating casing variation of a naturally fractured shale gas horizontal wellbore, which mainly includes the following steps:
[0054] 1) Establish a stress model of the formation-natural fracture-casing system before fracturing construction;
[0055] During the fracturing process of shale reservoirs, fracturing fluid may enter natural fractures due to fracture expansion or cement voids, causing the pressure in the fractures to rise. Once the pressure reaches a critical value, the sum of the fracture surface friction and the wellbore shear force is less than the formation slip shear force, which will cause natural fracture slip and further cause casing deformation. First, a force model of the formation-natural fracture-casing system before fracturing construction is established. From the plane stress analysis, it can be seen that the normal stress of the force unit in the x direction and the shear stress in the y direction before fracturing construction can be expressed as:
[0056]
[0057] Where: x , τ xy is the normal stress and shear stress of natural fractures, MPa; σ H , σ h are the maximum and minimum horizontal principal stresses, MPa; θ is the angle between the natural fracture and the maximum horizontal principal stress, rad.
[0058] 2) Establish the stress model for the case where the crack is not fully opened and the stress model for the case where the crack is fully opened;
[0059] (1) The crack is not fully expanded.
[0060] During fracturing, the fracturing fluid enters and supports the natural fracture. When the fracture is not fully opened, the main forces on the fracture surface in the x direction are the rock contact pressure and the fracturing fluid pressure, and the main forces on the fracture surface in the y direction are the friction force on the fracture surface and the shear force of the wellbore. Assuming that the force unit is stationary, a force balance relationship can be established in the x direction to obtain the shear stress of the wellbore to resist fracture slip under the current conditions:
[0061] τ c =(τ xy A f -f max ) / A c (2)
[0062] f max =μσ n A f (3)
[0063] Where: τ c is the shear stress on the casing, MPa; τ xy is the shear stress on the natural fracture, MPa; A f is the natural fracture area, m 2 ;f max is the friction force, N; A c is the cross-sectional area of the intersection of casing and natural fracture, m 2 ; μ is the friction coefficient of the crack surface; σ n is the rock contact normal stress on the fracture surface, MPa.
[0064] (2) The stress model of the crack being fully opened.
[0065] When the natural fracture is fully opened by the fracturing fluid, the rock masses on both sides of the fracture are no longer in contact, and the fracture surface is only acted on by pf in the x direction, and f is zero in the y direction. The shear force in the y direction is borne by the wellbore. Therefore, the wellbore shear stress under the current conditions can be calculated as:
[0066]
[0067] The cross-sectional area under the complete cement sheath is:
[0068]
[0069] The cross-sectional area of cement sheath failure is:
[0070]
[0071] Where: τ c is the shear stress on the casing, MPa; τ xy is the shear stress on the natural fracture, MPa; A f is the natural fracture area, m 2 ; A c is the cross-sectional area of the intersection of casing and natural fracture, m 2 ; R is the wellbore radius, m; α is the horizontal wellbore approach angle, °; r 1 is the inner diameter of the casing, m; r is the outer diameter of the casing, m; θ is the angle between the natural fracture and the maximum horizontal principal stress, rad.
[0072] 3) Considering the stress field changes caused by formation fluid injection, the corresponding pressure increment model is established;
[0073] The original in-situ stress is affected by the injection of fluid in the fractures of the previous fracturing section, resulting in stress redistribution. At this time, part of the in-situ stress is the change in the three-dimensional in-situ stress caused by the increase in pore pressure due to the injection of fluid in the fractures of the previous fracturing section, and the other part is the induced stress increment caused by the hydraulic fracture.
[0074] Correspondingly, the pore pressure increment model is established as;
[0075]
[0076] The critical pore pressure increment model for natural fracture activation and instability is established as follows:
[0077]
[0078] p(x,y)≥p 临界 (9)
[0079] Where: p(x,y) is the gap pressure increment at the plane (x,y); p fi,j (x, y) is the pore pressure generated by the jth segment of the i-th perforation cluster at the plane (x, y); p e is the original formation pressure, MPa; σ H , σ h are the maximum and minimum principal stresses, MPa.
[0080] The stress increment model of the change of three-dimensional ground stress due to the increase of injection pore pressure is:
[0081]
[0082] Δσ xy =σ 诱导mxy (12)
[0083] Where: ν is the Poisson's ratio of the reservoir rock, dimensionless; α is the Biot poroelastic coefficient, dimensionless; p p is the current formation pressure, MPa; σ 诱导mx , σ 诱导my , σ 诱导mxy are the normal stress and shear stress components generated by the crack in the x and y coordinates, MPa.
[0084] Because it is assumed that the casing variable is approximately equal to the tangential slip of the fracture, when establishing the casing slip DDM model, the stress field changes caused by fracturing are taken into account (i.e., the left side of equations (10, 11, 12) in this application): the original stress is affected by the injection of fluid in the fracture of the previous fracturing stage, resulting in stress redistribution. At this time, part of the ground stress is the increase in pore pressure in the fracture injection of the previous fracturing stage, which leads to changes in the three-dimensional ground stress (i.e., the first term on the right side of equations (10, 11, 12)); on the other hand, it is the induced stress increment generated by the hydraulic fracture (i.e., the second term on the right side of equations (10, 11, 12)). The pore pressure increment must be greater than the critical pore pressure increment for natural fracture activation instability before fracture tangential slip occurs.
[0085] It should be noted that in the original formation state, the fracture is in a relatively stable state, and the friction of the cross section is greater than the shear stress. As the fracturing fluid is continuously injected into the formation, the pore pressure continues to increase, and the friction of the cross section continues to decrease. When the friction is equal to the shear stress of the cross section, the fracture is in a critical state of sliding. With further injection of fracturing fluid, the fracture begins to lose stability and dislocate, acting on the casing, causing casing deformation.
[0086] 4) Use the DDM method to simulate the shear slip length of the natural fracture zone and calculate the set variable;
[0087] The corresponding model is:
[0088]
[0089] Where: is the shear stress and normal stress in the i-th fracture unit, MPa; p f is the fluid pressure in the natural fracture zone, MPa; β i is the angle between unit i and the positive x-axis, in rad; is the shear stress and normal stress generated by hydraulic fracture, MPa; is the shear stress caused by fracturing fluid loss, MPa.
[0090] Example 1
[0091] This example takes a well as an example to carry out the analysis of the main control impact of casing change. The basic parameters of the well are shown in Table 1.
[0092] Table 1 Basic parameters
[0093]
[0094] Based on the horizontal well fracturing sleeve variable calculation model established by the present invention, the sleeve variable main control influencing factors such as the natural fracture zone area, natural fracture zone friction coefficient, approach angle, natural fracture zone fluid pressure, shear stress, Young's modulus, Poisson's ratio and so on can be analyzed.
[0095] Based on the basic parameters of the well, the simulation of different natural fracture areas (10-100m 2 )The shear stress in the wellbore. Figure 2 The fracture is completely opened by the fracturing fluid, the natural fracture approach angle θ = 60°, and the fracture area (A) under different wellbore approach angle α values f ) and the wellbore shear stress. It can be seen that A f The larger the value, the c The larger the value, the c The value increases first and then decreases as the value of α increases. xy , A c When the value is constant, τ c Value A f The larger the natural fracture area, the higher the risk of wellbore shear failure.
[0096] Based on the basic parameters of the well, the wellbore approach angle is 80° and the natural fracture area is 10m 2 , the wellbore shear stress under the original formation pressure of 14.28MPa was simulated with friction coefficients of 0.6, 0.8, and 1 in different natural fracture zones. Figure 3 The crack is not fully opened by the fracturing fluid, α = 80°, p f =21MPa, A f =10m 2 When the crack friction coefficient (μ) is different, θ and τ c It can be seen that when the ground stress and the fracture geometric parameters are constant, the μ value increases, the maximum static friction of the fracture surface increases, and τ c Since the range of μ is limited, the change of μ has less influence on shear stress than p f The impact is significant.
[0097] Based on the basic parameters of the well, the friction coefficient of the natural fracture zone is 0.6, and the original formation pressure is 14.28MPa, the wellbore shear stress under different wellbore approach angles (0°-180°) and natural fracture zone approach angles (0°-180°) is simulated. Figure 4 It can be seen that when the crack direction is consistent with the maximum horizontal principal stress (σ H ) direction is parallel or perpendicular, the calculated value of wellbore shear stress is zero, which is independent of the wellbore approach angle α. The integrity of the cement sheath has a significant effect on the wellbore shear stress. The maximum wellbore shear stress under the condition of cement sheath failure (α=80°, θ=60°) is 8.75 times that under the condition of cement sheath integrity. The calculated value of wellbore shear stress is much higher than the shear strength of cement sheath (S s=11.5MPa), improving casing strength (S=178MPa) or cementing quality has limited effect on reducing the risk of casing deformation.
[0098] Based on the basic parameters of the well, combined with the natural fracture zone friction coefficient of 0.6, the wellbore approach angle of 80°, the natural fracture zone approach angle of 60°, the formation Young's modulus of 24900MPa, the Poisson's ratio of 0.19, and the original formation pressure of 14.28MPa, the wellbore shear stress under different fracture fluid pressures was simulated. Figure 5 It can be seen that when the natural fractures are not fully opened by the fracturing fluid, σ x -τ xy / μ <p f <σ x , the fluid pressure in the crack (p f ) is larger, the greater the wellbore shear stress is; f ≤σ x -τ xy / μ(15MPa), wellbore shear stress is 0; p f ≥σ x (29.9MPa), the wellbore shear stress is the largest.
[0099] Based on the basic parameters of the well, combined with the natural fracture zone friction coefficient of 0.6, the wellbore approach angle of 80°, the natural fracture zone approach angle of 60°, the formation Young's modulus of 24900MPa, the Poisson's ratio of 0.19, and the original formation pressure of 14.28MPa, the wellbore casing variables (slippage) at different wellbore positions under different shear stresses are simulated. Figure 6 It can be seen that the casing displacement (slippage) of the wellbore is symmetrical about the length direction of the natural fracture, that is, the casing displacement is small on both sides and large in the middle, forming an elliptical casing displacement; secondly, as the shear stress on the natural fracture increases from 0MPa to 15MPa, the casing displacement gradually increases, and the casing displacement increases the most when the shear stress increases from 5MPa to 10MPa. Therefore, in order to reduce the risk of casing displacement, on the one hand, the horizontal wellbore can be drilled at the end point of the natural fracture zone; on the other hand, the shear stress on the natural fracture zone can be reduced (which can be achieved by reducing the displacement or liquid volume).
[0100] Based on the basic parameters of the well, combined with the natural fracture zone friction coefficient of 0.6, the wellbore approach angle of 80°, the natural fracture zone approach angle of 60°, the Poisson's ratio of 0.19, and the original formation pressure of 14.28MPa, the wellbore casing variables (slippage) at different wellbore positions under different formation Young's moduli (14900MPa, 24900MPa and 34900MPa) were simulated. Figure 7It can be seen that the casing displacement (slippage) of the wellbore is symmetrical about the length direction of the natural fracture, that is, the casing displacement is small on both sides and large in the middle, forming an elliptical casing displacement; secondly, as the Young's modulus of the formation changes from 14900MPa to 34900MPa, the casing displacement gradually decreases, among which the casing displacement decreases the most when the Young's modulus changes from 14900MPa to 24900MPa.
[0101] Based on the basic parameters of the well, combined with the natural fracture zone friction coefficient of 0.6, the wellbore approach angle of 80°, the natural fracture zone approach angle of 60°, the formation Young's modulus of 24900MPa, and the original formation pressure of 14.28MPa, the wellbore casing variables (slippage) at different wellbore positions under different Poisson's ratios (0.19, 0.25 and 0.28) were simulated. Figure 8 It can be seen that the casing variation (slippage) of the wellbore is symmetrical about the length direction of the natural fracture, that is, the casing variation is small on both sides and large in the middle, which is an elliptical casing variation. Secondly, as the Poisson's ratio changes from 0.19 to 0.28, the casing variation gradually decreases. The change in Poisson's ratio does not affect the casing variation as much as the Young's modulus.
[0102] The main controlling factors of natural fracture area, friction coefficient of natural fracture zone, approach angle, fluid pressure in fracture, shear stress, Young's modulus and Poisson's ratio were analyzed by correlation analysis method. Fig. 9 From the formula, it can be seen that the order of the main controlling factors affecting the casing deformation is: fluid pressure in the fracture, Young's modulus, natural fracture area, approach angle, shear stress, friction coefficient of natural fracture zone, and Poisson's ratio.
[0103] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Claims
1. A method for calculating casing variables of horizontal wellbore for fracturing shale reservoirs with natural fractures, characterized in that: The following steps are involved: 1) Establish a stress model of the formation-natural fracture-casing system before fracturing construction; 2) Establish the stress model for the case where the crack is not fully opened and the stress model for the case where the crack is fully opened; 3) Considering the stress field changes caused by formation fracturing, the corresponding pressure increment model is established; 4) The DDM method is used to simulate the shear slip length of the natural fracture zone and perform set variable calculations.
2. The method for calculating casing variation of a horizontal wellbore according to claim 1, wherein the stress model of the formation-natural fracture-casing system in step 1) is: Where: x , τ xy is the normal stress and shear stress of natural fractures, MPa; σ H , σ h are the maximum and minimum horizontal principal stresses, MPa; θ is the angle between the natural fracture and the maximum horizontal principal stress, rad.
3. The method for calculating the horizontal wellbore casing variable according to claim 1, wherein the stress model of the fracture not fully propped open in step 2) is: t c =(t xy A f -f max ) / A c (2) f max =ms n A f (3) Where: τ c is the shear stress on the casing, MPa; τ xy is the shear stress on the natural fracture, MPa; A f is the natural fracture area, m 2 ;f max is the friction force, N; A c is the cross-sectional area of the intersection of casing and natural fracture, m 2 ; μ is the friction coefficient of the crack surface; σ n is the rock contact normal stress on the fracture surface, MPa.
4. The method for calculating the horizontal wellbore casing variable according to claim 1, wherein the stress model of the crack fully opened in step 2) is: The cross-sectional area under the complete cement sheath is: The cross-sectional area of cement sheath failure is: Where: τ c is the shear stress on the casing, MPa; τ xy is the shear stress on the natural fracture, MPa; A f is the natural fracture area, m 2 ; A c is the cross-sectional area of the intersection of casing and natural fracture, m 2 ; R is the wellbore radius, m; α is the horizontal wellbore approach angle, °; r1 is the inner diameter of the casing, m; r is the outer diameter of the casing, m; θ is the angle between the natural fracture and the maximum horizontal principal stress, rad.
5. The method for calculating horizontal wellbore casing variables according to claim 1, wherein step 3) further comprises: Establish pore pressure increment model; The critical pore pressure increment model for natural fracture activation and instability is: p(x,y)≥p 临界 (9) Where: p(x,y) is the gap pressure increment at plane (x,y); p fi,j (x, y) is the pore pressure generated by the jth segment of the i-th perforation cluster at the plane (x, y); p e is the original formation pressure, MPa; σ H , σ h are the maximum and minimum principal stresses, MPa.
6. The method for calculating the horizontal wellbore casing variable according to claim 1, wherein the step 3) further comprises: The stress increment model of the change of three-dimensional ground stress caused by the increase of injection pore pressure is: Board xy =s 诱导mxy (12) Where: ν is the Poisson's ratio of the reservoir rock, dimensionless; α is the Biot poroelastic coefficient, dimensionless; p p is the current formation pressure, MPa; σ 诱导mx , σ 诱导my , σ 诱导mxy are the normal stress and shear stress components generated by the crack in the x and y coordinates, MPa.
7. The method for calculating the casing variable of a horizontal wellbore according to claim 1, wherein the DDM method is used in step 4) to simulate the shear slip length of the natural fracture zone, and the corresponding model is: Where: is the shear stress and normal stress in the i-th fracture unit, MPa; p f is the fluid pressure in the natural fracture zone, MPa; β i is the angle between unit i and the positive x-axis, in rad; is the shear stress and normal stress generated by hydraulic fracture, MPa; is the shear stress caused by fracturing fluid loss, MPa.
8. A device for calculating casing variables of horizontal wellbore used for fracturing shale reservoirs with natural fractures, characterized in that: include: at least one processor and a memory in communication with the processor; The memory stores instructions that can be executed by the processor, and the instructions are executed by at least one processor to enable the at least one processor to execute the computing method described in any one of claims 1 to 7.
9. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the computing method according to any one of claims 1 to 7 are implemented.