A method for calculating a rock joint tip stress intensity factor

Abstract

The present application relates to a kind of rock crack tip stress intensity factor calculation method, belong to reservoir rock hydraulic fracturing field, comprising: obtaining the length l of conventional hydraulic fracturing crack and the pulse cycle hydraulic fracturing used fracturing fluid pump injection form, pulse frequency and amplitude;According to the length l of the conventional hydraulic fracturing crack and the pulse cycle hydraulic fracturing used fracturing fluid pump injection form, pulse frequency and amplitude, by means of transient flow model, the pressure wave velocity of fluid and the pressure distribution form in conventional hydraulic fracturing crack are obtained;According to the length l of the conventional hydraulic fracturing crack and the pressure wave velocity of fluid, the first order or second order inherent frequency of column in crack is calculated;According to the first order or second order inherent frequency of column in crack, the pulse cycle hydraulic fracturing used fracturing fluid pump injection form, pulse frequency and the pressure distribution form in the conventional hydraulic fracturing crack, obtain the rock crack tip stress intensity factor of three working conditions.

Description

Technical Field

[0001] This invention relates to the field of hydraulic fracturing of reservoir rocks, and specifically to a method for calculating the stress intensity factor at the rock fracture tip. Background Technology

[0002] Due to the low permeability of unconventional reservoirs such as low-permeability and tight reservoirs, hydraulic fracturing is necessary for their effective development to improve the flow of oil and gas within the reservoir. Compared with conventional hydraulic fracturing technology, pulsed circulation hydraulic fracturing technology can achieve more complex post-fracturing fracture morphologies with lower fracturing pressure. Both laboratory experiments and field tests have demonstrated the superiority and advanced nature of pulsed circulation hydraulic fracturing technology.

[0003] The calculation of the fracture tip stress intensity factor is a key parameter in fracture propagation simulation and reservoir rock fracturing pressure calculation. The fracture tip stress intensity factor can be obtained through analytical methods and numerical simulation. However, for calculating the fracture tip stress intensity factor in engineering applications, numerical simulation methods are cumbersome and difficult to converge due to the heterogeneity of the rock and the special morphology of the pressure distribution. Currently, there is no known analytical method for calculating the fracture tip stress intensity factor under unevenly distributed pressure.

[0004] In summary, current pulsed hydraulic fracturing methods lack analytical calculation methods for the fracture tip stress intensity factor, resulting in inaccurate design and optimization of pulsed circulating hydraulic fracturing construction parameters and poor pulsed circulating hydraulic fracturing effects. Summary of the Invention

[0005] To address the aforementioned problems, the purpose of this invention is to provide a method for calculating the rock fracture tip stress intensity factor. This method solves the problem that current pulsed hydraulic fracturing methods lack analytical calculation methods for the fracture tip stress intensity factor, leading to inaccurate design and optimization of pulsed circulating hydraulic fracturing construction parameters and poor pulsed circulating hydraulic fracturing results.

[0006] To achieve the above objectives, the present invention adopts the following technical solution:

[0007] This invention discloses a method for calculating the stress intensity factor at the tip of a rock fracture, comprising:

[0008] Step A: Obtain the length l of the conventional hydraulic fracturing fracture and the fracturing fluid injection method, pulse frequency, and amplitude used in pulsed circulation hydraulic fracturing. The fracturing fluid injection method used in pulsed circulation hydraulic fracturing includes pulse pressure and pulse flow rate.

[0009] Step B: Based on the length l of the conventional hydraulic fracturing fracture and the pumping method, pulse frequency and amplitude of the fracturing fluid used in the pulsed circulating hydraulic fracturing, the pressure wave velocity a of the fluid and the pressure distribution pattern P(x) in the conventional hydraulic fracturing fracture are obtained by using a transient flow model.

[0010] Step C: Calculate the first or second natural frequency of the fluid column within the fracture based on the length l of the conventional hydraulic fracturing fracture and the pressure wave velocity a of the fluid.

[0011] Step D: Based on the first or second natural frequency of the fluid column within the fracture, the pumping method of the fracturing fluid used in pulsed circulating hydraulic fracturing, the pulse frequency, and the pressure distribution pattern P(x) in the conventional hydraulic fracturing fracture, obtain the rock fracture tip stress intensity factor for three working conditions. Specifically, when the pumping method of the fracturing fluid used in pulsed circulating hydraulic fracturing is pulse pressure, if the input frequency is less than or equal to the first resonant frequency of the fluid column within the fracture, there is a peak amplitude located at the fracture tip, forming the first working condition; if the input frequency is less than or equal to the second resonant frequency of the fluid column within the fracture, there are two peak amplitudes located at one-third of the fracture length and at the fracture tip, forming the second working condition; when the pumping method of the fracturing fluid used in pulsed circulating hydraulic fracturing is pulse flow rate, the pressure distributed within the fracture is a uniform pressure equal at all locations, forming the third working condition.

[0012] Specifically, step B includes the following steps:

[0013] Step B1: Determine the formula for calculating the pressure wave velocity 'a' of the fluid. Its expression is:

[0014]

[0015] In the formula, K represents the elastic modulus of the liquid;

[0016] ρ is the density of the liquid;

[0017] E represents the elastic modulus of rock.

[0018] D represents the equivalent diameter of the crack.

[0019] e represents the roughness of the crack;

[0020] Step B2: Determine the formula for calculating the pressure distribution pattern P(x) in conventional hydraulic fracturing fractures. The expression is as follows:

[0021] P(x)=ρgH(x)

[0022] In the formula, g is the gravitational constant;

[0023] ρ is the density of the liquid;

[0024] H(x) is the average head of the cross section at position x of the fluid flow.

[0025] More specifically, step B2 includes the following steps:

[0026] Step B21: Based on the transient flow model, establish two partial differential equations that include the functional relationship between the cross-sectional average head H and the cross-sectional average velocity V;

[0027] Step B22: Transform the two partial differential equations established based on the transient flow model, including the functional relationship between the cross-sectional average head H and the cross-sectional average velocity V, into two sets of ordinary differential equations, namely the first characteristic line C. + The system of equations and the second characteristic line C - A system of equations, in which two systems of ordinary differential equations are obtained on characteristic lines with positive slopes and negative slopes, respectively. The system of equations on the characteristic lines with positive slopes is called the first characteristic line C. + A system of equations, with a negative slope, is called the second characteristic line C. - System of equations;

[0028] Step B23: For the first feature line C + The system of equations and the second characteristic line C - The equation system is meshed, and the flow Q of the left node A is established on the spatiotemporal graph after meshing through the left node A and the right node B. A The flow Q of right node B A And the pressure P at the next time node R R The associated iterative computation model;

[0029] Step B24: Solve for the pressure P at the next time node R. R The iterative calculation model sets boundary conditions;

[0030] Step B25: Based on the flow Q of left node A A The flow Q of right node B A And the pressure P at the next time node R R The associated iterative calculation model and its boundary conditions, combined with the length l of conventional hydraulic fracturing fractures, yield the pressure distribution pattern P(x) in conventional hydraulic fracturing fractures.

[0031] Specifically, the expressions for the two partial differential equations established based on the transient flow model, including the functional relationship between the cross-sectional average head H and the cross-sectional average velocity V, are as follows:

[0032]

[0033] In the formula, H is the cross-sectional average head H(x);

[0034] V is the average flow velocity across the cross section;

[0035] g is the gravitational constant;

[0036] x represents the location of the fluid flow;

[0037] t is the time of fluid flow;

[0038] τ0 is the shear stress coefficient;

[0039] ρ is the density of the liquid;

[0040] D is the equivalent diameter of the crack;

[0041] a represents the pressure wave velocity of the fluid.

[0042] Specifically, the first feature line C + The system of equations is expressed as follows:

[0043]

[0044] In the formula, V is the average flow velocity across the cross section;

[0045] t is the time of fluid flow;

[0046] g is the gravitational constant;

[0047] a is the pressure wave velocity of the fluid;

[0048] H is the average head of the cross section;

[0049] f is the Darcy-Weisbach friction coefficient;

[0050] D is the equivalent diameter of the crack;

[0051] x represents the location of the fluid flow;

[0052] |V| is the absolute value of the average flow velocity across the cross section;

[0053] The second characteristic line C - The system of equations is expressed as follows:

[0054]

[0055] In the formula, V is the average flow velocity across the cross section;

[0056] t is the time of fluid flow;

[0057] g is the gravitational constant;

[0058] a is the pressure wave velocity of the fluid;

[0059] H is the average head of the cross section;

[0060] f is the Darcy-Weisbach friction coefficient;

[0061] D is the equivalent diameter of the crack;

[0062] |V| is the absolute value of the average flow velocity across the cross section;

[0063] x represents the location of the fluid flow.

[0064] Specifically, step B23 includes the following steps:

[0065] Step B231: For the first feature line C + The system of equations and the second characteristic line C - The equations are meshed, and the cross-sectional average velocity V is transformed into flow rate Q through left node A and right node B on the meshed spatiotemporal graph. The flow rate Q of left node A is then obtained. A The flow Q of the right node B B ;

[0066] Step B232: Perform tests on the first feature line C respectively. + The system of equations and the second characteristic line C - The equations are subjected to difference processing to obtain the first characteristic line C. + The flow Q of the left node A A Pressure P at the next time node R R The equation expression and the second characteristic line C - The flow Q of the right node B A Pressure P at the next time node R R The equation expression together constitutes the flow Q of the left node A. A The flow Q of right node B A And the pressure P at the next time node R R The associated iterative computation model.

[0067] Specifically, the flow Q of the left node A A The flow Q of the right node B B The expressions are as follows:

[0068]

[0069] In the formula, A d The cross-sectional area of ​​the crack;

[0070] V A Let be the average flow velocity across the cross section of left node A;

[0071] V B Let be the average flow velocity at the cross section of right node B;

[0072] The flow Q of the left node A A The flow Q of right node B A And the pressure P at the next time node R R The expression for the associated iterative computation model is:

[0073]

[0074] In the formula, a is the pressure wave velocity of the fluid;

[0075] g is the gravitational constant;

[0076] A d The cross-sectional area of ​​the crack;

[0077] Q R For the traffic at the next time node R;

[0078] Q A Let A be the flow of the left node A;

[0079] P R The pressure on R at the next time point;

[0080] P A The pressure at the left node A;

[0081] f is the Darcy-Weisbach friction coefficient;

[0082] Δx is a variable representing the fluid position;

[0083] D is the equivalent diameter of the crack;

[0084] |Q A | represents the flow of the left node A;

[0085] Q B The flow of the right node B;

[0086] P B The pressure at the right node B;

[0087] |Q B | represents the flow of the right node B.

[0088] Specifically, step B24 includes the following specific steps:

[0089] The spatial domain start point is set at the wellhead location, i.e., x = 0, and the time domain start point is set at t = 0. The pumping method, pulse frequency, and amplitude of the fracturing fluid used in the pulsed circulating hydraulic fracturing are substituted into the flow rate Q of the left node A. A The flow Q of right node B A And the pressure P at the next time node R R In the associated iterative calculation model, the flow rate Q at the wellhead location as a function of time is obtained. x＝0 and pressure P x＝0 And the initial flow rate Q at all locations at time t=0. t＝0 and pressure P t＝0 .

[0090] Specifically, the expression for the first or second natural frequency of the liquid column within the slit is:

[0091]

[0092] In the formula, f n f1 represents the first or second natural frequency, n = 1 or 2, where f1 represents the first natural frequency and f2 represents the second natural frequency.

[0093] l represents the length of a conventional hydraulic fracturing fracture;

[0094] a represents the pressure wave velocity of the fluid.

[0095] Preferably, step D includes the following specific steps:

[0096] The relationship between the fracturing fluid injection method, pulse pressure frequency and first or second natural frequency used in the fracturing scheme is determined. The three working conditions are discussed, and the calculation formula of rock fracture tip stress intensity factor under each working condition is obtained.

[0097] First operating condition: When the fracturing fluid pumping method used in pulsed circulating hydraulic fracturing is pulsed pressure mode and the input frequency is less than or equal to the first natural frequency, the calculation formula for the fracture tip stress intensity factor under the first operating condition is obtained according to the pressure distribution pattern P(x) in the conventional hydraulic fracturing fracture:

[0098]

[0099] In the formula, K1 is the stress intensity factor at the seam tip under the first working condition;

[0100] l is the length of a conventional hydraulic fracturing fracture.

[0101] m represents the slope of the linear fit, and n represents the intercept of the linear fit. m and n are obtained by linearly fitting the pressure distribution pattern P(x) in conventional hydraulic fracturing fractures to P(x) = mx + n.

[0102] The second operating condition: When the fracturing fluid pumping method used in pulsed circulating hydraulic fracturing is pulsed pressure and the input frequency is less than or equal to the second natural frequency, the calculation formula for the fracture tip stress intensity factor under the second operating condition is obtained based on the pressure distribution pattern P(x) in the conventional hydraulic fracturing fracture.

[0103]

[0104] In the formula, K2 is the stress intensity factor at the seam tip under the second working condition;

[0105] l represents the length of a conventional hydraulic fracturing fracture;

[0106] n1 is the intercept parameter of the first segment after linear fitting;

[0107] n2 is the intercept parameter of the second segment after linear fitting;

[0108] n3 is the intercept parameter of the third segment after linear fitting;

[0109] The third operating condition: When the fracturing fluid pumping method used in pulsed circulating hydraulic fracturing is pulse flow rate and the pressure distribution within the fracture is uniformly distributed at all locations, based on the pressure distribution pattern P(x) in conventional hydraulic fracturing fractures, the calculation formula for the fracture tip stress intensity factor under the third operating condition is obtained:

[0110]

[0111] In the formula, K3 is the stress intensity factor at the seam tip under the third working condition;

[0112] l represents the length of a conventional hydraulic fracturing fracture;

[0113] U is a pressure constant, which is obtained by linearly fitting the pressure distribution pattern P(x) in a conventional hydraulic fracturing fracture to P(x) = U.

[0114] Compared with the prior art, the beneficial effects of the present invention are:

[0115] (I) This invention discloses a method for calculating the stress intensity factor at the rock fracture tip, including obtaining the length l of a conventional hydraulic fracturing fracture and the pumping method, pulse frequency, and amplitude of the fracturing fluid used in pulsed circulating hydraulic fracturing; obtaining the pressure wave velocity a of the fluid and the pressure distribution pattern P(x) in the conventional hydraulic fracturing fracture using a transient flow model; calculating the first or second natural frequency of the fluid column within the fracture; and obtaining the rock fracture tip stress intensity factor for three working conditions based on the first or second natural frequency of the fluid column within the fracture, the pumping method and pulse frequency of the fracturing fluid used in pulsed circulating hydraulic fracturing, and the pressure distribution pattern P(x) in the conventional hydraulic fracturing fracture.

[0116] (II) The rock fracture tip stress intensity factor disclosed in this invention, based on the different pressure distribution patterns formed in the hydraulic fracture by different fracturing fluid pumping methods in pulsed circulating hydraulic fracturing, provides expressions for the rock fracture tip stress intensity factor under three working conditions. This method can be used in pulsed circulating hydraulic fracturing to easily and accurately calculate the value of the rock fracture tip stress intensity factor, enhancing the accuracy of pulsed circulating hydraulic fracturing construction parameter design and optimization, improving the pulsed circulating hydraulic fracturing effect, and solving the problem that current pulsed hydraulic fracturing methods lack analytical calculation methods for the fracture tip stress intensity factor, leading to inaccurate pulsed circulating hydraulic fracturing construction parameter design and optimization, and poor pulsed circulating hydraulic fracturing effect. Attached Figure Description

[0117] Figure 1This is the spatiotemporal diagram after meshing provided in Example 1, using left node A and right node B to solve for the pressure P at the next time node R. R For example, where the horizontal axis represents the spatial domain and the vertical axis represents the time domain, i represents a node in the spatial domain, and j represents a node in the time domain, the pressure P at the next time node R is calculated. R During the process, both left node A and right node B can be considered as known nodes, and the next time node R is the unknown node to be solved. In the spatial domain, the next time node R is located to the right above left node A and to the left above right node B. In the time domain, left node A and right node B are at the same time step (e.g., t=0), and the next time node R is located at the next time step after left node A and right node B.

[0118] Figure 2 This is a diagram showing the relationship between the intra-crack pressure amplitude and the fracture location at the first and second natural frequencies provided in Embodiment 1 of the present invention. In the case of pumping fracturing fluid in the form of pulsed pressure, if the input frequency is less than or equal to the first resonant frequency, the obtained intra-crack pressure amplitude is the intra-crack pressure amplitude at the first natural frequency; if the input frequency is less than or equal to the second resonant frequency, the obtained intra-crack pressure amplitude is the intra-crack pressure amplitude at the second natural frequency.

[0119] Figure 3 This is a graph showing the relationship between the pressure amplitude inside the fracture and the fracture location when the input pulse pressure frequency is 11.05Hz under the condition of pumping fracturing fluid in the form of pulse pressure, as provided in Embodiment 1 of the present invention. Detailed Implementation

[0120] Exemplary embodiments of the invention will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the invention are shown in the drawings, it should be understood that the invention can be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided to enable a more thorough understanding of the invention and to fully convey the scope of the invention to those skilled in the art.

[0121] To address the lack of analytical calculation methods for fracture tip stress intensity factor in current pulsed hydraulic fracturing (PHF) methods, which leads to inaccurate design and optimization of PHF construction parameters and poor PHF fracturing performance, this invention discloses a method for calculating the rock fracture tip stress intensity factor. Based on the different pressure distribution patterns formed in the hydraulic fracture by different fracturing fluid injection methods in PHF, expressions for the rock fracture tip stress intensity factor under three operating conditions are provided. This method can be used in PHF to easily and accurately calculate the value of the rock fracture tip stress intensity factor, enhancing the accuracy of PHF construction parameter design and optimization, improving PHF fracturing performance, and solving the problem of inaccurate design and optimization of PHF construction parameters and poor PHF fracturing performance caused by the lack of analytical calculation methods for fracture tip stress intensity factor in current PHF methods.

[0122] Example 1: A method for calculating the stress intensity factor at the tip of a rock fracture.

[0123] Example 1 provides a method for calculating the stress intensity factor at the tip of a rock fracture, comprising the following steps:

[0124] Step A: Obtain the length l of the conventional hydraulic fracturing fracture and the pumping method, pulse frequency, and amplitude of the fracturing fluid used in pulsed circulation hydraulic fracturing;

[0125] Pulsed circulation hydraulic fracturing is often combined with conventional hydraulic fracturing. After obtaining simple, long, straight hydraulic fractures through conventional hydraulic fracturing, the length of the conventional hydraulic fracturing fractures is obtained using methods such as microseismic monitoring.

[0126] The design of the pulsed circulating hydraulic fracturing scheme yields the fracturing fluid injection method, pulse frequency, and amplitude. The fracturing fluid injection method for pulsed circulating hydraulic fracturing includes pulse pressure and pulse flow rate.

[0127] Step B: Based on the length of the conventional hydraulic fracturing fracture and the pumping method, pulse frequency, and amplitude of the fracturing fluid used in the pulsed circulating hydraulic fracturing, the pressure wave velocity a and the pressure distribution pattern P(x) in the conventional hydraulic fracturing fracture are obtained using a transient flow model. This specifically includes the following steps:

[0128] Step B1: Determine the formula for calculating the pressure wave velocity 'a' of the fluid. Its expression is:

[0129]

[0130] In the formula, K represents the elastic modulus of the liquid;

[0131] ρ is the density of the liquid;

[0132] E represents the elastic modulus of rock.

[0133] D represents the equivalent diameter of the crack.

[0134] e represents the roughness of the crack.

[0135] Step B2: Determine the formula for calculating the pressure distribution pattern P(x) in conventional hydraulic fracturing fractures. The expression is as follows:

[0136] P(x)=ρgH(x) (Equation 2)

[0137] In the formula, g is the gravitational constant;

[0138] ρ is the density of the liquid;

[0139] H(x) is the average head of the cross section at position x of the fluid flow, which is a pressure value.

[0140] As a specific implementation method, the cross-sectional average head H(x) is obtained using a transient flow model, including the following steps:

[0141] Step B21: Based on the transient flow model, establish two partial differential equations that include the functional relationship between the cross-sectional average head H and the cross-sectional average velocity V. Their expressions are as follows:

[0142]

[0143] In the formula, H is the cross-sectional average head H(x);

[0144] V is the average flow velocity across the cross section;

[0145] g is the gravitational constant;

[0146] x represents the location of the fluid flow;

[0147] t is the time of fluid flow;

[0148] τ0 is the shear stress coefficient;

[0149] ρ is the density of the liquid;

[0150] D is the equivalent diameter of the crack;

[0151] a represents the pressure wave velocity of the fluid.

[0152] Specifically, the expression for τ0, the shear stress coefficient, is as follows:

[0153]

[0154] In the formula, ρ is the liquid density;

[0155] V is the average flow velocity across the cross section;

[0156] |V| is the absolute value of the average flow velocity across the cross section;

[0157] f is the Darcy-Weisbach friction coefficient.

[0158] Specifically, the expression for the Darcy-Weisbach friction coefficient f is:

[0159]

[0160] In the formula, e represents the roughness of the crack;

[0161] D represents the equivalent diameter of the crack;

[0162] Re represents the Reynolds coefficient.

[0163] Specifically, the expression for the Reynolds coefficient Re is:

[0164]

[0165] In the formula, ρ is the liquid density;

[0166] V is the average flow velocity across the cross section;

[0167] D is the equivalent diameter of the crack;

[0168] μ is the viscosity of the liquid.

[0169] Specifically, the expression for the equivalent diameter D of the crack is:

[0170] D = 4 × (A d / S) (Formula 7)

[0171] In the formula, A d The cross-sectional area of ​​the crack;

[0172] S is the perimeter of the cross section.

[0173] In order to obtain the cross-sectional average head H and cross-sectional average velocity V of the transient flow model, the two partial differential equations of the transient flow model need to be processed.

[0174] Step B22: Transform the two partial differential equations established based on the transient flow model, including the functional relationship between the cross-sectional average head H and the cross-sectional average velocity V, into two sets of ordinary differential equations, namely the first characteristic line C. + The system of equations (as shown in Equation 8) and the second characteristic line C - The system of equations (as shown in Equation 9) consists of two sets of ordinary differential equations obtained on characteristic lines with positive and negative slopes, respectively. The system of equations on the characteristic lines with positive slopes is called the first characteristic line C. +A system of equations, with a negative slope, is called the second characteristic line C. - A system of equations.

[0175] Wherein, the first feature line C + The system of equations is expressed as follows:

[0176]

[0177] In the formula, V is the average flow velocity across the cross section;

[0178] t is the time of fluid flow;

[0179] g is the gravitational constant;

[0180] a is the pressure wave velocity of the fluid;

[0181] H is the average head of the cross section;

[0182] f is the Darcy-Weisbach friction coefficient;

[0183] D is the equivalent diameter of the crack;

[0184] x represents the location of the fluid flow;

[0185] |V| represents the absolute value of the average flow velocity across the cross section.

[0186] Wherein, the second characteristic line C - The system of equations is expressed as follows:

[0187]

[0188] In the formula, V is the average flow velocity across the cross section;

[0189] t is the time of fluid flow;

[0190] g is the gravitational constant;

[0191] a is the pressure wave velocity of the fluid;

[0192] H is the average head of the cross section;

[0193] f is the Darcy-Weisbach friction coefficient;

[0194] D is the equivalent diameter of the crack;

[0195] |V| is the absolute value of the average flow velocity across the cross section;

[0196] x represents the location of the fluid flow.

[0197] Step B23: For the first feature line C + The system of equations and the second characteristic line C -The equation system is meshed, and the pressure P for solving the next time node R is established on the spatiotemporal diagram after meshing through left node A and right node B. R The iterative computation model includes the following steps:

[0198] Step B231: For the first feature line C + The system of equations and the second characteristic line C - The equations are meshed, and on the resulting spatiotemporal graph, the cross-sectional average velocity V is transformed into a flow rate Q through left node A and right node B. The flow rate Q at left node A is then used to represent the flow rate V. A The flow Q of the right node B B The expressions are as follows:

[0199]

[0200] In the formula, A d The cross-sectional area of ​​the crack;

[0201] V A Let be the average flow velocity across the cross section of left node A;

[0202] V B Let be the average flow velocity across the cross section of right node B.

[0203] Step B232: Perform tests on the first feature line C respectively. + The system of equations and the second characteristic line C - The equations are subjected to difference processing to obtain the first characteristic line C. + The flow Q of the left node A A Pressure P at the next time node R R The equation expression and the second characteristic line C - The flow Q of the right node B A Pressure P at the next time node R R The equation expression together constitutes the flow Q of the left node A. A The flow Q of right node B A And the pressure P at the next time node R R The associated iterative computation model, its expression:

[0204]

[0205] In the formula, a is the pressure wave velocity of the fluid;

[0206] g is the gravitational constant;

[0207] A d The cross-sectional area of ​​the crack;

[0208] Q RFor the traffic at the next time node R;

[0209] Q A Let A be the flow of the left node A;

[0210] P R The pressure on R at the next time point;

[0211] P A The pressure at the left node A;

[0212] f is the Darcy-Weisbach friction coefficient;

[0213] Δx is a variable representing the fluid position;

[0214] D is the equivalent diameter of the crack;

[0215] |Q A | represents the flow of the left node A;

[0216] Q B The flow of the right node B;

[0217] P B The pressure at the right node B;

[0218] |Q B | represents the flow of the right node B.

[0219] Step B24: Solve for the pressure P at the next time node R. R The iterative calculation model sets boundary conditions, specifically including the following steps:

[0220] refer to Figure 1 The spatial domain start point is set at the wellhead location, i.e., x = 0, and the time domain start point is set at t = 0. The pumping method, pulse frequency, and amplitude of the fracturing fluid used in the pulsed circulating hydraulic fracturing are substituted into the flow rate Q of the left node A. A The flow Q of right node B A And the pressure P at the next time node R R In the associated iterative calculation model, the flow rate Q at the wellhead location as a function of time is obtained. x＝0 and pressure P x＝0 And the initial flow rate Q at all locations at time t=0. t＝0 and pressure P t＝0 .

[0221] For example, when the fracturing fluid pumping method used in the pulsed circulating hydraulic fracturing is pulsed pressure pumping, the pressure P at the wellhead position changes with time. x＝0 =20sin(10πt), then the pressure P at the left node A at this time is... A=20sin(10πt). When the fracturing fluid pumping method used in the pulsed circulating hydraulic fracturing is pulse flow rate pumping, the flow rate Q at the wellhead position changes with time. x＝0 =2sin(4πt), then the flow rate Q at the left node A is at this time. A =2sin(4πt).

[0222] Step B25: Based on the flow Q of left node A A The flow Q of right node B A And the pressure P at the next time node R R The associated iterative calculation model and its boundary conditions, combined with the length l of conventional hydraulic fracturing fractures, yield the pressure distribution pattern P(x) in conventional hydraulic fracturing fractures.

[0223] Specifically, starting from x = 0, ΔZ, 2ΔZ, 3ΔZ, ... are continuously calculated at fixed intervals ΔZ until x equals the length l of the conventional hydraulic fracturing fracture. The pressure at each location on the conventional hydraulic fracturing fracture can be obtained. The pressure at each location acts perpendicularly on the conventional hydraulic fracturing fracture; therefore, the pressure at each location on the conventional hydraulic fracturing fracture represents the pressure distribution pattern P(x) within the conventional hydraulic fracturing fracture.

[0224] Studies have found that when the fracturing fluid pumping method used in pulsed circulating hydraulic fracturing is pulse pressure, the pressure distribution within the fracture is a non-uniform amplitude pulse pressure. This non-uniform pressure distribution is due to the formation of standing waves, the shape of which is influenced by the input pulse pressure, amplitude, and fracture length. When the pressure within a fracture of known length is non-uniformly distributed along the fracture length, the number of peak pulse pressure amplitudes is determined by the relationship between the input frequency and the resonant frequency of the fluid column within the fracture. If the input frequency is less than or equal to the first-order resonant frequency of the fluid column within the fracture, there is one peak amplitude located at the fracture tip, forming the first operating condition. If the input frequency is less than or equal to the second-order resonant frequency of the fluid column within the fracture, there are two peak amplitudes located at one-third of the fracture length and at the fracture tip, forming the second operating condition. Figure 2 As shown.

[0225] When the fracturing fluid pump used in pulsed circulating hydraulic fracturing is in pulse flow form, the pressure distributed in the fracture is a uniform pressure equal at all locations, forming the third working condition.

[0226] Step C: Calculate the first or second natural frequency of the fluid column within the fracture based on the length of the conventional hydraulic fracturing fracture and the pressure wave velocity a of the fluid.

[0227] The expression for the first or second natural frequency of the liquid column within the slit is as follows:

[0228]

[0229] In the formula f n f1 represents the first or second natural frequency, n = 1 or 2, where f1 represents the first natural frequency and f2 represents the second natural frequency.

[0230] l represents the length of a conventional hydraulic fracturing fracture;

[0231] a represents the pressure wave velocity of the fluid.

[0232] Step D: Based on the first or second natural frequency of the fluid column within the fracture, the pumping method of the fracturing fluid used in pulsed circulating hydraulic fracturing, the pulse frequency, and the pressure distribution pattern P(x) in the conventional hydraulic fracturing fracture, obtain the rock fracture tip stress intensity factor for three operating conditions. Specifically, when the pumping method of the fracturing fluid used in pulsed circulating hydraulic fracturing is pulse pressure, if the input frequency is less than or equal to the first resonant frequency of the fluid column within the fracture, there is a peak amplitude located at the fracture tip, forming the first operating condition; if the input frequency is less than or equal to the second resonant frequency of the fluid column within the fracture, there are two peak amplitudes located at one-third of the fracture length and at the fracture tip, forming the second operating condition; when the pumping method of the fracturing fluid used in pulsed circulating hydraulic fracturing is pulse flow rate, the pressure distributed within the fracture is a uniform pressure equal at all locations, forming the third operating condition, including the following specific steps:

[0233] The relationship between the fracturing fluid injection method, pulse pressure frequency and first or second natural frequency used in the fracturing scheme is determined. The three working conditions are discussed, and the calculation formula of rock fracture tip stress intensity factor under each working condition is obtained.

[0234] First operating condition: When the fracturing fluid pumping method used in pulsed circulating hydraulic fracturing is pulsed pressure mode and the input frequency is less than or equal to the first natural frequency, the calculation formula for the fracture tip stress intensity factor under the first operating condition is obtained according to the pressure distribution pattern P(x) in the conventional hydraulic fracturing fracture:

[0235]

[0236] In the formula, K1 is the stress intensity factor at the seam tip under the first working condition;

[0237] l is the length of a conventional hydraulic fracturing fracture.

[0238] m represents the slope of the linear fit, and n represents the intercept of the linear fit. m and n are obtained by linearly fitting the pressure distribution pattern P(x) in conventional hydraulic fracturing fractures to P(x) = mx + n.

[0239] The second operating condition: When the fracturing fluid pumping method used in pulsed circulating hydraulic fracturing is pulsed pressure and the input frequency is less than or equal to the second natural frequency, the calculation formula for the fracture tip stress intensity factor under the second operating condition is obtained based on the pressure distribution pattern P(x) in the conventional hydraulic fracturing fracture.

[0240]

[0241] In the formula, K2 is the stress intensity factor at the seam tip under the second working condition;

[0242] l represents the length of a conventional hydraulic fracturing fracture;

[0243] n1 is the intercept parameter of the first segment after linear fitting;

[0244] n2 is the intercept parameter of the second segment after linear fitting;

[0245] n3 is the intercept parameter of the third segment after linear fitting.

[0246] Under this condition, the pressure distribution pattern P(x) in a conventional hydraulic fracturing fracture can be linearly fitted in three segments as follows: P(x) = m1x + n1, P(x) = m2x + n2, and P(x) = m3x + n3. Here, m1 is the slope parameter of the first segment after linear fitting; m2 is the slope parameter of the second segment after linear fitting; and m3 is the slope parameter of the third segment after linear fitting.

[0247] It is worth mentioning that the frequency used in controlling fracturing fluid injection with pulse pressure is generally less than or equal to the second-order natural frequency of the fluid column, so the formula for calculating the stress intensity factor at the fracture tip at a higher frequency is not considered.

[0248] For example, refer to Figure 3 Taking a 100m half-fracture hydraulic fracture subjected to pulsed circulation hydraulic fracturing as an example, fracturing fluid is pumped in a pulsed pressure manner with a pulse frequency of 11.05Hz. The pressure amplitude fractions at various locations within the fracture are obtained using a transient flow model. Figure 3 As shown. Piecewise linear fitting yields:

[0249] P = 0.3708x + 1.4

[0250] P = -0.3531x + 24.237

[0251] P = 0.3235x - 9.858

[0252] The first and second natural frequencies of the crack are calculated to be 6.6 Hz and 11.05 Hz, respectively. This falls under the second working condition. Based on the stress intensity factor at the crack tip under the second working condition, the stress intensity factor can be calculated to be 1.8369.

[0253] The third operating condition: When the fracturing fluid pumping method used in pulsed circulating hydraulic fracturing is pulse flow rate and the pressure distribution within the fracture is uniformly distributed at all locations, based on the pressure distribution pattern P(x) in conventional hydraulic fracturing fractures, the calculation formula for the fracture tip stress intensity factor under the third operating condition is obtained:

[0254]

[0255] In the formula, K3 is the stress intensity factor at the seam tip under the third working condition;

[0256] l represents the length of a conventional hydraulic fracturing fracture;

[0257] U is a pressure constant, which is obtained by linearly fitting the pressure distribution pattern P(x) in a conventional hydraulic fracturing fracture to P(x) = U.

[0258] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for calculating the stress intensity factor at the tip of a rock fracture, characterized in that, include: Step A: Obtain the length of a conventional hydraulic fracturing fracture. The pumping method, pulse frequency, and amplitude of the fracturing fluid used in pulsed circulating hydraulic fracturing, wherein the pumping method of the fracturing fluid used in pulsed circulating hydraulic fracturing includes pulse pressure form and pulse flow rate form; Step B: Based on the length of the conventional hydraulic fracturing fracture... The pressure wave velocity of the fluid is obtained using a transient flow model, taking into account the pumping method, pulse frequency, and amplitude of the fracturing fluid used in the pulsed circulating hydraulic fracturing. a Pressure distribution patterns in conventional hydraulic fracturing fractures ; Step C: Based on the length of the conventional hydraulic fracturing fracture... and the pressure wave velocity of the fluid a The first or second natural frequency of the liquid column within the slit is calculated using the following expression: In the formula, It is a first-order or second-order natural frequency. =1 or 2, where Represents the first-order natural frequency. Indicates the second-order natural frequency; This refers to the length of a conventional hydraulic fracturing fracture. a The pressure wave velocity of the fluid; Step D: Based on the first or second natural frequency of the fluid column within the fracture, the pumping method of the fracturing fluid used in the pulsed circulating hydraulic fracturing, the pulse frequency, and the pressure distribution pattern in the conventional hydraulic fracturing fracture. The stress intensity factor at the rock fracture tip was obtained for three operating conditions. Specifically, when the fracturing fluid pumping method used in pulsed circulating hydraulic fracturing is pulse pressure mode, if the input frequency is less than or equal to the first-order resonant frequency of the fluid column within the fracture, there is a peak amplitude located at the fracture tip, forming the first operating condition. If the input frequency is less than or equal to the second-order resonant frequency of the fluid column within the fracture, there are two peak amplitudes located at one-third of the fracture length and at the fracture tip, forming the second operating condition. When the fracturing fluid pumping method used in pulsed circulating hydraulic fracturing is pulse flow mode, the pressure distributed within the fracture is a uniform pressure equal at all locations, forming the third operating condition. The specific steps are as follows: The relationship between the fracturing fluid injection method, pulse pressure frequency and first or second natural frequency used in the fracturing scheme is determined. The three working conditions are discussed, and the calculation formula of rock fracture tip stress intensity factor under each working condition is obtained. The first operating condition: When the fracturing fluid pump used in pulsed circulating hydraulic fracturing is in pulsed pressure mode and the input frequency is less than or equal to the first natural frequency, according to the pressure distribution pattern in the conventional hydraulic fracturing fracture. The formula for calculating the stress intensity factor at the seam tip under the first working condition is obtained as follows: In the formula, This is the stress intensity factor at the seam tip under the first working condition; l The length of a conventional hydraulic fracturing fracture m This represents the slope of the linear fit. n Represents the linear fit intercept. m and n Pressure distribution patterns in conventional hydraulic fracturing fractures Linear fitting =mx+n get; The second operating condition: When the fracturing fluid pump used in pulsed circulating hydraulic fracturing is in pulsed pressure mode and the input frequency is less than or equal to the second-order natural frequency, based on the pressure distribution pattern in the conventional hydraulic fracturing fracture... The formula for calculating the stress intensity factor at the seam tip under the second working condition is obtained as follows: In the formula, This refers to the stress intensity factor at the seam tip under the second working condition. l This refers to the length of a conventional hydraulic fracturing fracture. The intercept parameter for the first segment after linear fitting; The intercept parameter for the second segment after linear fitting; The intercept parameter for the third segment after linear fitting; The third operating condition: When the fracturing fluid pump used in pulsed circulating hydraulic fracturing is in pulse flow mode and the pressure distribution within the fracture is uniformly distributed at all locations, based on the pressure distribution pattern in conventional hydraulic fracturing fractures. The formula for calculating the stress intensity factor at the seam tip under the third working condition is obtained as follows: In the formula, This refers to the stress intensity factor at the seam tip under the third working condition. l This refers to the length of a conventional hydraulic fracturing fracture. U The pressure constant is determined by the pressure distribution pattern in conventional hydraulic fracturing fractures. After linear fitting, =U .

2. The method for calculating the rock fracture tip stress intensity factor according to claim 1, characterized in that, Step B includes the following steps: Step B1: Determine the pressure wave velocity of the fluid. a The formula for calculating it is expressed as follows: In the formula, K Indicates the elastic modulus of a liquid; The density of the liquid; E Indicates the elastic modulus of rock. D Indicates the equivalent diameter of the crack. e Indicates the roughness of the crack; Step B2: Determine the pressure distribution pattern in conventional hydraulic fracturing fractures The formula for calculating it is expressed as follows: In the formula, It is the gravitational constant; The density of the liquid; Location of fluid flow The average head at the cross section.

3. The method for calculating the stress intensity factor at the rock fracture tip according to claim 2, characterized in that, Step B2 includes the following steps: Step B21: Based on the transient flow model, establish a cross-sectional average head. H and cross-sectional average velocity Two partial differential equations relating functions; Step B22: Based on the transient flow model, the cross-sectional average head is calculated... H and cross-sectional average velocity The two partial differential equations relating functions are transformed into two systems of ordinary differential equations, which are the first characteristic lines. System of equations and second characteristic line A system of equations, in which two systems of ordinary differential equations are obtained on characteristic lines with positive slopes and negative slopes, respectively. The system of equations on the characteristic lines with positive slopes is called the first characteristic line. A system of equations, with a negative slope, is called the second characteristic line. System of equations; Step B23: For the first feature line The system of equations and the second characteristic line The equation system is meshed, and the flow of the left node A is established on the spatiotemporal graph after meshing through the left node A and the right node B. Flow of right node B and the pressure at the next time node R The associated iterative computation model; Step B24: Solve for the pressure at the next time node R. The iterative calculation model sets boundary conditions; Step B25: Based on the flow of left node A Flow of right node B and the pressure at the next time node R The associated iterative computational model and its boundary conditions, as well as the length of conventional hydraulic fracturing fractures. l Obtain the pressure distribution morphology in conventional hydraulic fracturing fractures. .

4. The method for calculating the rock fracture tip stress intensity factor according to claim 3, characterized in that, The transient flow model was established, including the cross-sectional average head. H and cross-sectional average velocity The expressions for the two partial differential equations relating the functions are: In the formula, H The average head of the cross section ; The average flow velocity across the cross section; It is the gravitational constant; The location of fluid flow; t The time of fluid flow; It is the shear stress coefficient; ρ The density of the liquid; The equivalent diameter of the crack; a The pressure wave velocity of the fluid.

5. The method for calculating the rock fracture tip stress intensity factor according to claim 3, characterized in that, First feature line The system of equations is expressed as follows: In the formula, The average flow velocity across the cross section; t The time of fluid flow; It is the gravitational constant; a The pressure wave velocity of the fluid; H The average head of the cross section; The Darcy-Weisbach friction coefficient; The equivalent diameter of the crack; x The location of fluid flow; This represents the absolute value of the average flow velocity across the cross section. Second feature line The system of equations is expressed as follows: In the formula, The average flow velocity across the cross section; t The time of fluid flow; It is the gravitational constant; a The pressure wave velocity of the fluid; H The average head of the cross section; The Darcy-Weisbach friction coefficient; The equivalent diameter of the crack; This represents the absolute value of the average flow velocity across the cross section. x This indicates the location of the fluid flow.

6. The method for calculating the stress intensity factor at the rock fracture tip according to claim 5, characterized in that, Step B23 includes the following steps: Step B231: For the first feature line The system of equations and the second characteristic line The equations are meshed, and the average flow velocity of the cross section is calculated on the spatiotemporal graph through left node A and right node B. Deformation is performed to increase the average flow velocity of the cross-section. Transform into flow Q Obtain the flow of the left node A respectively. and the flow of right node B ; Step B232: Perform tests on the first feature lines respectively. The system of equations and the second characteristic line The equations are subjected to difference processing to obtain the first characteristic lines respectively. The flow of the left node A Pressure at the next time node R The equation expression and the second characteristic line The flow of the right node B Pressure at the next time node R The equation expression together constitutes the flow of the left node A. Flow of right node B and the pressure at the next time node R The associated iterative computation model.

7. The method for calculating the rock fracture tip stress intensity factor according to claim 6, characterized in that, The flow of the left node A and the flow of right node B The expressions are as follows: In the formula, A d The cross-sectional area of ​​the crack; V A Let be the average flow velocity at the cross section of left node A; V B Let be the average flow velocity at the cross section of right node B; The flow of the left node A Flow of right node B and the pressure at the next time node R The expression for the associated iterative computation model is: In the formula, a The pressure wave velocity of the fluid; It is the gravitational constant; The cross-sectional area of ​​the crack; For the next time node R Traffic; Left node A Traffic; For the next time node R The pressure; Left node A The pressure; The Darcy-Weisbach friction coefficient; The variable represents the fluid position; For cracks equivalent diameter; Left node A Traffic; Right node B Traffic; Right node B The pressure; Let B be the flow of the right node.

8. The method for calculating the rock fracture tip stress intensity factor according to claim 7, characterized in that, Step B24 includes the following specific steps: The starting point of the spatial domain is set at the wellhead location, i.e. x At point =0, the time domain start point is set to time t=0. The pumping method, pulse frequency, and amplitude of the fracturing fluid used in the pulsed circulating hydraulic fracturing are substituted into the flow rate at left node A. Flow of right node B and the pressure at the next time node R In the associated iterative calculation model, the flow rate at the wellhead location as a function of time is obtained. Q x=0 and pressure P x=0 and the initial flow rate at all locations at time t=0. Q t=0 and pressure P t=0 .

Images