A non-isodensity flow-limiting perforation optimization design method

Abstract

The invention discloses a non-isodensity flow-limiting perforation optimization design method, which belongs to the technical field of hydraulic fracturing, and comprises the following steps: step 1, calculating the minimum horizontal geostress of a target fracturing section according to well logging data; step 2, calculating the induced stress along the direction of the minimum horizontal geostress generated by adjacent fractured sections on the target fracturing section; step 3, considering the mutual interference between hydraulic fractures in the target fracturing section, calculating the induced stress generated by the expansion of hydraulic fractures in the target fracturing section; step 4, calculating the geostress magnitude of each perforation cluster in the target fracturing section; step 5, formulating several initial perforation schemes according to the geostress magnitude relationship calculated in step 4; step 6, calculating the flow distribution results under different initial perforation schemes, and determining the optimal perforation scheme according to the principle of minimizing the difference in flow distribution between clusters; the method of the invention provides a more scientific and accurate optimization technical means for the perforation parameter design of staged multi-cluster fracturing of horizontal wells.

Description

Technical Field

[0001] The invention belongs to the technical field of hydraulic fracturing, and in particular relates to a non-isodensity flow-limiting perforation optimization design method. Background Art

[0002] As an efficient measure to increase the production of oil and gas wells, the multi-cluster fracturing technology of horizontal wells has been widely used in the development of unconventional oil and gas reservoirs. This technology forms multiple fracture clusters by implementing fracturing on multiple segments of horizontal wells, thereby effectively increasing the contact area between oil and gas and the wellbore and improving the efficiency of oil and gas flow. As a key step in the fracturing process, the main purpose of perforation is to create a channel between the wellbore and the oil and gas layer to provide a path for fracturing fluid injection and oil and gas flow. However, in the actual fracturing process, multi-cluster fracturing of horizontal wells often faces the problem of difficult balanced development of multiple clusters of fractures, which limits the improvement of oil and gas well production capacity.

[0003] In order to promote the balanced development of multiple clusters of fractures, the flow-limited perforation technology was proposed. Flow-limited perforation increases the friction resistance of the perforation holes by reducing the hole diameter and the number of holes, and increases the bottom hole pressure, forcing the fracturing fluid to be evenly distributed among the perforation clusters. However, the existing isodensity flow-limited perforation method may lead to excessive friction resistance of the holes and excessive construction pressure by uniformly reducing the number or diameter of holes in each perforation cluster, greatly increasing the difficulty and cost of construction, and thus affecting the economic benefits of oil and gas wells.

[0004] In order to overcome this limitation, non-isodensity flow-limiting perforation technology was introduced. Non-isodensity flow-limiting perforation technology optimizes the flow distribution of perforation clusters based on the heterogeneity of ground stress and the interaction between fractures. However, the current non-isodensity flow-limiting perforation schemes are based on engineering experience and trial and error, lacking systematic theoretical support and accurate calculation methods, which limits the optimization accuracy and applicability of non-isodensity flow-limiting perforation technology, thereby reducing the fracturing effect. Therefore, there is an urgent need for a non-isodensity flow-limiting perforation optimization design method to improve the accuracy and applicability of the perforation scheme. Summary of the invention

[0005] The purpose of the present invention is to provide a non-isodensity flow-limiting perforation optimization design method, which analyzes the dynamic distribution of flow between clusters by comprehensively considering factors such as ground stress heterogeneity, inter-segment stress interference and intra-segment stress interference, and realizes fine optimization of the perforation scheme.

[0006] To achieve the above object, the present invention provides a non-isodensity flow-limiting perforation optimization design method, comprising the following steps:

[0007] Step 1: Calculate the minimum horizontal ground stress of the target fracturing section based on the well logging data;

[0008] Step 2: Considering the influence of existing hydraulic fractures on newly developed hydraulic fractures, calculate the induced stress along the horizontal minimum ground stress direction generated by the adjacent fractured section on the target fractured section;

[0009] Step 3: Considering the mutual interference between the hydraulic fractures in the target fracturing section, calculate the induced stress generated by the expansion of the hydraulic fractures in the target fracturing section;

[0010] Step 4: Calculate the ground stress of each perforation cluster in the target fracturing section;

[0011] Step 5: Formulate several initial perforation plans based on the ground stress relationship calculated in step 4;

[0012] Step 6: Calculate the flow distribution results under different initial perforation schemes, and determine the optimal perforation scheme based on the principle of minimizing the difference in flow distribution between clusters.

[0013] Preferably, in step 1, the specific calculation expression of the minimum horizontal in-situ stress of the target fracturing section is as follows:

[0014]

[0015] In the formula, σ Hmin is the horizontal minimum ground stress, in MPa; ξ1 is the horizontal maximum stress structural coefficient, dimensionless; ξ2 is the horizontal minimum stress structural coefficient, dimensionless; v is the static Poisson's ratio, dimensionless; E is the static Young's modulus, in MPa; α is the Biot coefficient, dimensionless; p s represents the crack induced stress value;

[0016] The determination of static Young's modulus E and static Poisson's ratio v is based on the following steps: first, experimental measurements are conducted on core samples from adjacent wells to obtain the baseline values ​​of static rock mechanics parameters; second, the corresponding dynamic rock mechanics parameters are obtained through well logging technology; finally, the obtained parameters are used to accurately fit the required static Young's modulus E and static Poisson's ratio v using the regression analysis method; the specific calculation formulas for static Young's modulus E and static Poisson's ratio v are as follows:

[0017] E=a×E d +b (2)

[0018] ν=c×μ d +d (3)

[0019] Where a, b, c, and d are all regression coefficients, which are related to the regional geological background and are dimensionless; E d is the dynamic Young's modulus, in MPa; μ d is the dynamic Poisson's ratio, dimensionless;

[0020] Among them, the dynamic Young's modulus of rock Ed and the dynamic Poisson's ratio μ d It is calculated based on the P-wave, S-wave time difference and density logging data. The specific expression is as follows:

[0021]

[0022] In the formula, Δt s is the shear wave time difference, in μs / m; Δt p is the longitudinal wave time difference, in μs / m; ρ b is the rock density in g / cm 3 ;

[0023] The specific calculation formula of Biot coefficient α is as follows:

[0024]

[0025] In the formula, ρ ma is the maximum rock density, in g / cm 3 ; V p , V s are the longitudinal wave velocity and the transverse wave velocity, both in m / μs; V map , V mas is the maximum longitudinal wave velocity and transverse wave velocity, both in m / μs.

[0026] Preferably, the specific calculation formula for the induced stress generated by the adjacent fractured sections on the target fractured section along the horizontal minimum geostress direction in step 2 is as follows:

[0027]

[0028] In the formula, σ ss m is the induced stress at any unit point i, in MPa; A ss m,j A is the plane strain elastic coefficient parallel to the crack wall, dimensionless; sn m,j is the plane strain elastic coefficient perpendicular to the crack wall, dimensionless; D s j is the tangential displacement discontinuity of the jth micro-segment of the fracture unit, in m; D n j is the normal displacement discontinuity on the jth micro-segment of the crack unit, in m; G m,j is the fracture height correction factor, dimensionless; M is the total number of units into which the hydraulic fracture is divided; where the fracture height correction factor G m,j The calculation formula is as follows:

[0029]

[0030] Where h is the height of the hydraulic fracture, in meters; d m,j is the distance from crack unit m to crack unit j, in m.

[0031] Preferably, the calculation process of the induced stress generated by the expansion of the hydraulic fracture in the target fracturing section in step 3 is as follows:

[0032] First, it is roughly assumed that the flow distribution of each fracture is uniform, and the net pressure p of the fracture is solved based on the PKN model:

[0033]

[0034] Secondly, considering that the stress interference between cracks decays with the increase of crack spacing, the coefficient g of the induced stress value decaying with distance is:

[0035]

[0036] Then the crack induced stress value p is estimated s for:

[0037] p s =cgp (11)

[0038] Where E is the static Young's modulus, in MPa; μ is the viscosity, in MPa·s; Q is the injection flow rate, in m 3 / min, N is the number of perforation clusters in a single stage, dimensionless; v is Poisson's ratio, dimensionless; h r is the height of hydraulic fracture, in m, d c is the cluster spacing, in meters; c is the empirical position correction coefficient, where c=0.5 for the outermost perforation cluster position in the segment and c=1 for other perforation cluster positions.

[0039] Preferably, the calculation expression of the ground stress of each perforation cluster in the target fracturing stage in step 4 is as follows:

[0040] p sum =σ Hmin +σ ss m +p s (12)

[0041] Preferably, the specific process of formulating several initial perforation schemes according to the geostress magnitude relationship calculated in step 4 in step 5 is as follows:

[0042] S51, selecting a cluster with representativeness and clear geological characteristics from each perforation cluster in the fracturing stage as a reference cluster, and assigning an initial perforation number based on the geological characteristics and ground stress level of the cluster;

[0043] S52. For the remaining perforation clusters, the friction coefficient and the number of perforations are calculated according to the stress difference between them and the reference cluster to achieve non-isodensity perforation. The specific expression is as follows:

[0044]

[0045] In the formula, λ i is the friction coefficient of the i-th perforation cluster, dimensionless; ΔF is the stress difference, in MPa; Q is the injection flow rate, in m 3 / min; N is the number of perforation clusters in a single stage; λ0 is the friction coefficient of the reference cluster, dimensionless; n i is the number of perforations in the ith perforation cluster, which is an integer; ρ is the density of the fracturing fluid, in kg / m 3 ;d i is the perforation hole diameter of the ith perforation cluster, in m; K di is the perforation discharge coefficient of the ith perforation cluster, dimensionless;

[0046] S53, changing the initial perforation number of the reference cluster, recalculating the perforation number of the remaining clusters according to the difference in ground stress, and formulating several different perforation schemes.

[0047] Preferably, in step 6, the flow distribution results under different initial perforation schemes are calculated, and the specific process of determining the optimal perforation scheme according to the principle of minimizing the difference in flow distribution between clusters is as follows:

[0048] S61, initialize traffic distribution, the specific formula is as follows:

[0049]

[0050] In the formula, Q i is the flow rate of each cluster, in m 3 / min; Q t is the total flow rate, in m 3 / min; N is the number of perforation clusters in a single stage;

[0051] S62, construct the Jacobian matrix J and the residual function F to update the flow distribution, where the Jacobian matrix is ​​an N×N square matrix, and the elements J on the diagonal ii represents the flow rate Q of each perforation cluster i Friction coefficient of the hole itself λ i The contribution of the off-diagonal elements J i,i+1 Represents the flow influence between adjacent perforation clusters. The specific representation of the Jacobian matrix is ​​as follows:

[0052]

[0053] The residual function is an N-dimensional vector representing the flow rate Q of each perforation cluster. i The difference between the expected value under given conditions is expressed as follows:

[0054]

[0055] Where F is the residual function, dimensionless; λ i is the friction coefficient of the ith perforation cluster, dimensionless; Q i is the flow rate of the ith perforation cluster, in m 3 / min; F i is the i-th residual function component, dimensionless; F n is the nth residual function component, dimensionless; σ i Hmin is the horizontal minimum principal stress of the i-th perforation cluster; Q t is the total flow rate, in m 3 / min; σ i+1 Hmin is the horizontal minimum principal stress of the i+1th perforation cluster; i+1 is the friction coefficient of the i+1th perforation cluster, dimensionless;

[0056] S63, by solving the linear equation group JΔQ=F to obtain the flow correction amount ΔQ, and update the current flow distribution. The specific expression is as follows:

[0057] ΔQ=J -1 F (18)

[0058] S64, continue to iterate until the following convergence conditions are met, and determine the optimal perforation scheme based on the principle of minimizing the difference in flow distribution between clusters. The specific expression is as follows:

[0059] |ΔQ|<10 -10 (19)

[0060] Therefore, the present invention adopts the above-mentioned non-isodensity flow-limiting perforation optimization design method, which can give the optimal perforation scheme by iteratively optimizing the flow distribution based on the comprehensive consideration of the heterogeneity of ground stress, the stress interference between segments and within segments. This method can not only improve the production capacity of oil and gas wells, but also reduce construction costs and enhance the economic benefits of oil and gas wells, which has important practical significance and application value for promoting the efficient development of oil and gas fields.

[0061] The technical solution of the present invention is further described in detail below through the accompanying drawings and embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] Figure 1 It is an overall flow chart of a non-isodensity flow-limiting perforation optimization design method of the present invention. DETAILED DESCRIPTION

[0063] The following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the invention claimed for protection, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0064] See also Figure 1 , a non-isodensity limiting perforation optimization design method, comprising the following steps:

[0065] Step 1: Calculate the minimum horizontal ground stress of the target fracturing section based on the well logging data; the specific calculation expression is as follows:

[0066]

[0067] In the formula, σ Hmin is the horizontal minimum ground stress, in MPa; ξ1 is the horizontal maximum stress structural coefficient, dimensionless; ξ2 is the horizontal minimum stress structural coefficient, dimensionless; v is the static Poisson's ratio, dimensionless; E is the static Young's modulus, in MPa; α is the Biot coefficient, dimensionless; p s represents the crack induced stress value;

[0068] The determination of static Young's modulus E and static Poisson's ratio v is based on the following steps: first, experimental measurements are conducted on core samples from adjacent wells to obtain the baseline values ​​of static rock mechanics parameters; second, the corresponding dynamic rock mechanics parameters are obtained through well logging technology; finally, the obtained parameters are used to accurately fit the required static Young's modulus E and static Poisson's ratio v using the regression analysis method; the specific calculation formulas for static Young's modulus E and static Poisson's ratio v are as follows:

[0069] E=a×E d +b (2)

[0070] ν=c×μ d +d (3)

[0071] Where a, b, c, and d are all regression coefficients, which are related to the regional geological background and are dimensionless; E d is the dynamic Young's modulus, in MPa; μ d is the dynamic Poisson's ratio, dimensionless;

[0072] Among them, the dynamic Young's modulus of rock E d and the dynamic Poisson's ratio μ d It is calculated based on the P-wave, S-wave time difference and density logging data. The specific expression is as follows:

[0073]

[0074] In the formula, Δt s is the shear wave time difference, in μs / m; Δt p is the longitudinal wave time difference, in μs / m; ρ b is the rock density in g / cm 3 ;

[0075] The specific calculation formula of Biot coefficient α is as follows:

[0076]

[0077] In the formula, ρ ma is the maximum rock density, in g / cm 3 ; V p , V s are the longitudinal wave velocity and the transverse wave velocity, both in m / μs; V map , V mas is the maximum longitudinal wave velocity and transverse wave velocity, both in m / μs.

[0078] Step 2: Considering the influence of existing hydraulic fractures on newly developed hydraulic fractures, calculate the induced stress in the direction of the minimum horizontal ground stress generated by the adjacent fractured section on the target fractured section; the specific calculation formula is as follows:

[0079]

[0080] In the formula, σ ss m is the induced stress at any unit point i, in MPa; A ss m,j A is the plane strain elastic coefficient parallel to the crack wall, dimensionless; sn m,j is the plane strain elastic coefficient perpendicular to the crack wall, dimensionless; D s j is the tangential displacement discontinuity of the jth micro-segment of the fracture unit, in m; D n j is the normal displacement discontinuity on the jth micro-segment of the crack unit, in m; G m,j is the fracture height correction factor, dimensionless; M is the total number of units into which the hydraulic fracture is divided; where the fracture height correction factor G m,j The calculation formula is as follows:

[0081]

[0082] Where h is the height of the hydraulic fracture, in meters; d m,j is the distance from crack unit m to crack unit j, in m.

[0083] Step 3: Considering the mutual interference between the hydraulic fractures in the target fracturing section, calculate the induced stress generated by the expansion of the hydraulic fractures in the target fracturing section; the specific calculation process is as follows:

[0084] First, it is roughly assumed that the flow distribution of each fracture is uniform, and the net pressure p of the fracture is solved based on the PKN model:

[0085]

[0086] Secondly, considering that the stress interference between cracks decays with the increase of crack spacing, the coefficient g of the induced stress value decaying with distance is:

[0087]

[0088] The crack induced stress value is estimated ps for:

[0089] p s =cgp (11)

[0090] Where, E is the static Young's modulus, in MPa; μ is the viscosity, in MPa·s; Q is the injection flow rate, in m 3 / min, N is the number of perforation clusters in a single stage, dimensionless; v is Poisson's ratio, dimensionless; h r is the height of hydraulic fracture, in m, d c is the cluster spacing, in meters; c is the empirical position correction coefficient, where c=0.5 for the outermost perforation cluster position in the segment and c=1 for other perforation cluster positions.

[0091] Step 4: Calculate the ground stress of each perforation cluster in the target fracturing section; the specific calculation expression is as follows:

[0092] p sum =σ Hmin +σ ss m +p s (12)

[0093] Step 5: Formulate several initial perforation plans based on the ground stress relationship calculated in step 4. The specific process is as follows:

[0094] S51, selecting a cluster with representativeness and clear geological characteristics from each perforation cluster in the fracturing stage as a reference cluster, and assigning an initial perforation number based on the geological characteristics and ground stress level of the cluster;

[0095] S52. For the remaining perforation clusters, the friction coefficient and the number of perforations are calculated according to the stress difference between them and the reference cluster to achieve non-isodensity perforation. The specific expression is as follows:

[0096]

[0097] In the formula, λ i is the friction coefficient of the i-th perforation cluster, dimensionless; ΔF is the stress difference, in MPa; Q is the injection flow rate, in m 3 / min; N is the number of perforation clusters in a single stage; λ0 is the friction coefficient of the reference cluster, dimensionless; n i is the number of perforations in the ith perforation cluster, which is an integer; ρ is the density of the fracturing fluid, in kg / m 3 ;d i is the perforation hole diameter of the ith perforation cluster, in m; K di is the perforation discharge coefficient of the ith perforation cluster, dimensionless;

[0098] S53, changing the initial perforation number of the reference cluster, recalculating the perforation number of the remaining clusters according to the difference in ground stress, and formulating several different perforation schemes.

[0099] Step 6: Calculate the flow distribution results under different initial perforation schemes, and determine the optimal perforation scheme based on the principle of minimizing the difference in flow distribution between clusters. The specific process is as follows:

[0100] S61, initialize traffic distribution, the specific formula is as follows:

[0101]

[0102] In the formula, Q i is the flow rate of each cluster, in m 3 / min; Q t is the total flow rate, in m 3 / min; N is the number of perforation clusters in a single stage;

[0103] S62, construct the Jacobian matrix J and the residual function F to update the flow distribution, where the Jacobian matrix is ​​an N×N square matrix, and the elements J on the diagonal ii represents the flow rate Q of each perforation cluster i Friction coefficient of the hole itself λ i The contribution of the off-diagonal elements J i,i+1 Represents the flow influence between adjacent perforation clusters. The specific representation of the Jacobian matrix is ​​as follows:

[0104]

[0105] The residual function is an N-dimensional vector representing the flow rate Q of each perforation cluster. i The difference between the expected value under given conditions is expressed as follows:

[0106]

[0107] Where F is the residual function, dimensionless; λ i is the friction coefficient of the ith perforation cluster, dimensionless; Q i is the flow rate of the ith perforation cluster, in m 3 / min; F i is the i-th residual function component, dimensionless; F n is the nth residual function component, dimensionless; σ i Hmin is the horizontal minimum principal stress of the i-th perforation cluster; Q t is the total flow rate, in m 3 / min; σ i+1 Hmin is the horizontal minimum principal stress of the i+1th perforation cluster; i+1 is the friction coefficient of the i+1th perforation cluster, dimensionless;

[0108] S63, by solving the linear equation group JΔQ=F to obtain the flow correction amount ΔQ, and update the current flow distribution. The specific expression is as follows:

[0109] ΔQ=J -1 F (18)

[0110] S64, continue to iterate until the following convergence conditions are met, and determine the optimal perforation scheme based on the principle of minimizing the difference in flow distribution between clusters. The specific expression is as follows:

[0111] |ΔQ|<10 -10 (19)

[0112] Example

[0113] Taking the 6th section of the shale oil horizontal well ZM6 in the eastern region as an example, the construction parameters of this section are shown in Table 1 below. Table 1 Construction parameters

[0114] Basic parameters Value Basic parameters Value <![CDATA[Injection flow rate (m 3 / min)]]> 16 Number of perforation clusters per stage 6 Fracturing fluid viscosity (mPa·s) 5 Hydraulic fracture height (m) 30 Young's modulus (MPa) 30000 <![CDATA[Fracturing fluid density (kg / m 3 )]]> 1050 Cluster spacing (m) 10 Hole flow coefficient 0.95 Construction time (min) 100 Poisson's ratio (dimensionless) 0.2 Initial perforation diameter (m) 0.0095

[0115] The specific process is as follows:

[0116] Step 1: According to the field logging data of Well ZM6, the minimum horizontal stress of each perforation cluster in the sixth section is calculated to be 49MPa, 49MPa, 52MPa, 51.5MPa, 50.5MPa, and 49MPa respectively;

[0117] Step 2, calculate that the induced stresses of the fractured section on each cluster along the horizontal minimum ground stress direction are 3.3 MPa, 2.8 MPa, 2.4 MPa, 2.1 MPa, 1.9 MPa, and 1.8 MPa respectively;

[0118] Step 3, calculate the stress interference values ​​within the segment as 1.12MPa, 2.24MPa, 2.24MPa, 2.24MPa, 2.24MPa, 1.12MPa;

[0119] Step 4: Calculate the in-situ stress of the perforation clusters in the fracturing section, and obtain the stress distribution of 53.42MPa, 54.04MPa, 56.64MPa, 55.84MPa, 54.64MPa, and 51.92MPa.

[0120] Step 5: According to the relationship between the in-situ stresses among the perforation clusters, set up a variety of feasible initial perforation schemes. Select cluster 1 as the reference cluster, and set the initial perforation numbers to 8, 9, and 10 based on its in-situ stress. Use the stress difference value to determine the number of perforations in other perforation clusters to obtain the following three initial perforation schemes.

[0121] Scheme 1: The initial perforation number of the reference cluster is 8, and the perforation numbers of the other clusters are calculated based on the stress difference and are 9, 11, 10, 9, and 8, respectively.

[0122] Solution 2: The initial perforation number of the reference cluster is adjusted to 9, and the perforation numbers of the other clusters are calculated to be 10, 11, 10, 11, and 9 according to the stress difference.

[0123] Solution 3: The initial perforation number of the reference cluster is adjusted to 10, and the perforation numbers of the other clusters are calculated to be 11, 12, 9, 11, and 10 according to the stress difference.

[0124] Step 6: Iteratively optimize the flow distribution of perforation. According to the principle of minimizing the difference in flow distribution between clusters, the optimal perforation scheme is finally determined to be Scheme 1, that is, the number of perforations in each perforation cluster in the 6th section of Well ZM6 is 9, 10, 11, 10, 11, and 9, respectively.

[0125] According to the non-isodensity flow-limiting perforation optimization design method provided by the present invention, the perforation scheme of ZM6 well was optimized. The daily production of ZM6 well after fracturing was 16.8×10 4 m 3 , which is 1.48 times the production of the adjacent well ZM7, indicating that the method of the present invention can effectively improve the ultimate recovery rate of oil and gas wells by accurately calculating and optimizing flow distribution.

[0126] Therefore, the present invention adopts the above-mentioned non-isodensity flow-limiting perforation optimization design method. Compared with the traditional perforation technology, the core advantage is that it can comprehensively consider the heterogeneity of ground stress, the stress interference between segments and within segments, and give the optimal perforation scheme by iteratively optimizing the flow distribution. This method not only improves the production capacity of oil and gas wells, but also reduces construction costs and enhances the economic benefits of oil and gas wells. It has important practical significance and application value for promoting the efficient development of oil and gas fields, and provides a more scientific and accurate optimization technology for the design of horizontal well segmented multi-cluster fracturing perforation parameters.

[0127] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solution of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solution to deviate from the spirit and scope of the technical solution of the present invention.

Claims

1. A non-isodensity flow-limiting perforation optimization design method, characterized in that: The following steps are involved: Step 1: Calculate the minimum horizontal ground stress of the target fracturing section based on the well logging data; Step 2: Considering the influence of existing hydraulic fractures on newly developed hydraulic fractures, calculate the induced stress along the horizontal minimum ground stress direction generated by the adjacent fractured section on the target fractured section; Step 3: Considering the mutual interference between the hydraulic fractures in the target fracturing section, calculate the induced stress generated by the expansion of the hydraulic fractures in the target fracturing section; Step 4: Calculate the ground stress of each perforation cluster in the target fracturing section; Step 5: Formulate several initial perforation plans based on the ground stress relationship calculated in step 4; Step 6: Calculate the flow distribution results under different initial perforation schemes, and determine the optimal perforation scheme based on the principle of minimizing the difference in flow distribution between clusters; In step 6, the flow distribution results under different initial perforation schemes are calculated. According to the principle of minimizing the difference in flow distribution between clusters, the specific process of determining the optimal perforation scheme is as follows: S61, initialize traffic distribution, the specific formula is as follows: In the formula, Q i is the flow rate of each cluster, in m 3 / min; Q t is the total flow rate, in m 3 / min; N is the number of perforation clusters in a single stage; S62, construct the Jacobian matrix J and the residual function F to update the flow distribution, where the Jacobian matrix is ​​an N×N square matrix, and the elements J on the diagonal ii represents the flow rate Q of each perforation cluster i Friction coefficient of the hole itself λ i The contribution of the off-diagonal elements J i,i+1 Represents the flow influence between adjacent perforation clusters. The specific representation of the Jacobian matrix is ​​as follows: The residual function is an N-dimensional vector representing the flow rate Q of each perforation cluster. i The difference between the expected value under given conditions is expressed as follows: Where F is the residual function, dimensionless; λ i is the friction coefficient of the ith perforation cluster, dimensionless; Q i is the flow rate of the ith perforation cluster, in m 3 / min; F i is the i-th residual function component, dimensionless; F n is the nth residual function component, dimensionless; σ i Hmin is the horizontal minimum principal stress of the i-th perforation cluster; Q t is the total flow rate, in m 3 / min; σ i+1 Hmin is the horizontal minimum principal stress of the i+1th perforation cluster; i+1 is the friction coefficient of the i+1th perforation cluster, dimensionless; S63, by solving the linear equation group JΔQ=F to obtain the flow correction amount ΔQ, and update the current flow distribution. The specific expression is as follows: ΔQ=J -1 F (18) S64, continue to iterate until the following convergence conditions are met, and determine the optimal perforation scheme based on the principle of minimizing the difference in flow distribution between clusters. The specific expression is as follows: |ΔQ|<10 -10 (19) 2. The non-isodensity flow-limiting perforation optimization design method according to claim 1 is characterized in that: In step 1, the specific calculation expression of the minimum horizontal ground stress of the target fracturing section is as follows: In the formula, σ Hmin is the horizontal minimum ground stress, in MPa; ξ1 is the horizontal maximum stress structural coefficient, dimensionless; ξ2 is the horizontal minimum stress structural coefficient, dimensionless; v is the static Poisson's ratio, dimensionless; E is the static Young's modulus, in MPa; α is the Biot coefficient, dimensionless; p s represents the crack induced stress value; The determination of static Young's modulus E and static Poisson's ratio v is based on the following steps: first, experimental measurements are conducted on core samples from adjacent wells to obtain the baseline values ​​of static rock mechanics parameters; second, the corresponding dynamic rock mechanics parameters are obtained through well logging technology; finally, the obtained parameters are used to accurately fit the required static Young's modulus E and static Poisson's ratio v using the regression analysis method; the specific calculation formulas for static Young's modulus E and static Poisson's ratio v are as follows: E=a×E d +b (2) n=c×μ d +d (3) Where a, b, c, and d are all regression coefficients, which are related to the regional geological background and are dimensionless; E d is the dynamic Young's modulus, in MPa; μ d is the dynamic Poisson's ratio, dimensionless; Among them, the dynamic Young's modulus of rock E d and the dynamic Poisson's ratio μ d It is calculated based on the P-wave, S-wave time difference and density logging data. The specific expression is as follows: In the formula, Δt s is the shear wave time difference, in μs / m; Δt p is the longitudinal wave time difference, in μs / m; ρ b is the rock density in g / cm 3 ; The specific calculation formula of Biot coefficient α is as follows: In the formula, ρ ma is the maximum rock density, in g / cm 3 ; V p , V s are the longitudinal wave velocity and the transverse wave velocity, both in m / μs; V map , V mas is the maximum longitudinal wave velocity and transverse wave velocity, both in m / μs.

3. The non-isodensity flow-limiting perforation optimization design method according to claim 2 is characterized in that: The specific calculation formula for the induced stress generated by the adjacent fractured section on the target fractured section along the horizontal minimum ground stress direction in step 2 is as follows: In the formula, σ ss m is the induced stress at any unit point i, in MPa; A ss m,j is the plane strain elastic coefficient parallel to the crack wall, dimensionless; A sn m,j is the plane strain elastic coefficient perpendicular to the crack wall, dimensionless; D s j is the tangential displacement discontinuity of the jth micro-segment of the fracture unit, in m; D n j is the normal displacement discontinuity on the jth micro-segment of the crack unit, in m; G m,j is the seam height correction factor, dimensionless; M is the total number of units into which the hydraulic fractures are divided; among them, the fracture height correction factor G m,j The calculation formula is as follows: Where h is the height of hydraulic fracture, in meters; d m,j is the distance from crack unit m to crack unit j, in m.

4. The non-isodensity flow-limiting perforation optimization design method according to claim 3 is characterized in that: The calculation process of the induced stress generated by the expansion of hydraulic fractures in the target fracturing section in step 3 is as follows: First, it is roughly assumed that the flow distribution of each fracture is uniform, and the net pressure p of the fracture is solved based on the PKN model: Secondly, considering that the stress interference between cracks decays with the increase of crack spacing, the coefficient g of the induced stress value decaying with distance is: Then the crack induced stress value p is estimated s for: p s =cgp (11) Where, E is the static Young's modulus, in MPa; μ is the viscosity, in MPa·s; Q is the injection flow rate, in m 3 / min, N is the number of perforation clusters in a single stage, dimensionless; v is Poisson's ratio, dimensionless; h r is the height of hydraulic fracture, in m, d c is the cluster spacing, in meters; c is the empirical position correction coefficient, where c=0.5 for the outermost perforation cluster position in the segment and c=1 for other perforation cluster positions.

5. The non-isodensity flow-limiting perforation optimization design method according to claim 4 is characterized in that: The calculation expression of the ground stress of each perforation cluster in the target fracturing section in step 4 is as follows: p sum =s Hmin +s ss m +p s (12)。 6. The non-isodensity flow-limiting perforation optimization design method according to claim 5 is characterized in that: The specific process of formulating several initial perforation schemes in step 5 according to the ground stress relationship calculated in step 4 is as follows: S51, selecting a cluster from each perforation cluster in the fracturing stage as a reference cluster, and assigning an initial perforation number based on the geological characteristics and ground stress level of the cluster; S52. For the remaining perforation clusters, the friction coefficient and the number of perforations are calculated according to the stress difference between them and the reference cluster. The specific expression is as follows: In the formula, λ i is the friction coefficient of the ith perforation cluster, dimensionless; ΔF is the stress difference, in MPa; Q is the injection flow rate, in m 3 / min; N is the number of perforation clusters in a single stage; λ0 is the friction coefficient of the reference cluster, dimensionless; n i is the number of perforations in the ith perforation cluster, which is an integer; ρ is the density of the fracturing fluid, in kg / m 3 ; d i is the perforation hole diameter of the ith perforation cluster, in m; K di is the perforation discharge coefficient of the ith perforation cluster, dimensionless; S53, changing the initial perforation number of the reference cluster, recalculating the perforation number of the remaining clusters according to the difference in ground stress, and formulating several different perforation schemes.

Images