A shale reservoir fracturing casing deformation calculation method

Abstract

This invention discloses a method for calculating casing deformation in shale reservoir fracturing, relating to the field of shale oil and gas development technology. The method includes: integrating microseismic, logging, core, and drilling data to construct a fracture network-geomechanical coupled model; calculating the direction cosine matrix using fracture triaxial stress and dip angle, assessing activation risk, and screening out high-risk fractures; solving for slip using triangular mesh boundary elements as displacement boundaries, iteratively calculating casing variables using an eight-node hexahedral elastoplastic casing mesh, and plotting contour maps; and iterating the model in real-time with newly added microseismic data during construction, dynamically optimizing perforation cluster spacing, fracturing scale, displacement, and well spacing. This invention fully utilizes a large amount of field microseismic data to achieve fine characterization of reservoir fractures and slip risk analysis, and completes multi-point casing variable prediction in the wellbore, providing key technical support for fracturing scheme decision-making and risk avoidance in deep shale gas development.

Description

Technical Field

[0001] This application relates to the field of shale oil and gas development technology, and in particular to a method for calculating the deformation of fracturing casing in shale reservoirs. Background Technology

[0002] During fracturing operations in deep shale gas formations, complex geological conditions and fracturing parameters lead to significant casing deformation (hereinafter referred to as casing deformation) and even hydraulic channeling. Casing deformation affects wellbore integrity, causing frequent instances of lost sections in horizontal wells, and even resulting in some wells losing production altogether. Hydraulic channeling is more challenging, causing reduced production in numerous adjacent wells, severely impacting the fracturing process and production ramp-up. Existing research indicates that the primary factor inducing casing deformation is fluid-induced fault shear slippage, which in turn causes casing deformation due to the shearing action on the casing. Based on the degree of casing inner diameter deformation, casing deformation is classified into five levels, and corresponding fracturing stimulation schemes can be implemented according to the deformation level to mitigate its impact. Therefore, accurate calculation and prediction of casing deformation are crucial for deep shale gas development.

[0003] Existing technologies mainly consider the fracture morphology as an idealized two-dimensional vertical fracture and use two-dimensional Mohr's circles to solve for activation risk and slip. However, the existing technologies predict the deformation of the fracturing casing in shale reservoirs in an overly idealized and inaccurate manner, thus failing to accurately reflect the casing deformation risk.

[0004] Therefore, there is an urgent need for a method to calculate the deformation of fracturing casing in shale reservoirs, so as to effectively improve the prediction accuracy of fracturing casing deformation in shale reservoirs. Summary of the Invention

[0005] Therefore, it is necessary to provide a method for calculating the deformation of fracturing casing in shale reservoirs to address the aforementioned technical problems.

[0006] The present invention adopts the following technical solution:

[0007] This invention provides a method for calculating the deformation of fracturing casing in shale reservoirs, including:

[0008] Acquire microseismic data, imaging logging data, core observation data, and drilling and logging data generated in shale reservoirs during fracturing operations; construct a fracture network model based on the microseismic data, imaging logging data, and core observation data;

[0009] Based on the fracture network model, well logging data, and drilling data, a geomechanical model is constructed; the geomechanical model includes the magnitude and direction of the three principal stresses on each fracture surface;

[0010] Calculate the direction cosine matrix based on the angle between the normal of each fracture surface and the direction of the three principal stresses; calculate the activation risk value of each fracture surface based on the magnitude of the three principal stresses and the elements in the direction cosine matrix, and designate the fracture surfaces with activation risk values ​​greater than a preset threshold as high-risk fracture surfaces.

[0011] The high-risk fracture surface is meshed into three-dimensional triangular surface mesh elements. A system of linear algebraic equations is established with the continuous invariant shear displacement of the three-dimensional triangular surface mesh elements as unknowns. The three-dimensional slip distribution of the high-risk fracture surface is obtained by solving the equations.

[0012] Using the three-dimensional slip distribution as the displacement boundary condition, a linear hardening elastoplastic constitutive equation is introduced, and an eight-node hexahedral element is used to construct a casing mesh model. By iteratively calculating the casing mesh model, the nonlinear deformation of the casing is obtained, and a casing inner diameter deformation cloud map is plotted.

[0013] Preferably, the method further includes: iteratively updating the fracture network model and the geomechanical model based on newly added microseismic data, and dynamically optimizing the parameters of perforation cluster spacing, fracturing scale / discharge and well spacing by combining the fracturing fracture length, stress disturbance range and casing deformation risk prediction results.

[0014] Preferably, a fracture network model is constructed based on microseismic data, imaging logging data, and core observation data, specifically including:

[0015] Based on microseismic data, the main fracture trend surface and attitude information of the horizontal section of the fractured well are determined; the attitude information includes: fracture surface dip angle, fracture surface dip direction, and fracture surface area.

[0016] An initial fracture network model for the work area was established based on microseismic data, imaging logging data, and core observation data.

[0017] Based on the imaging logging data and the core observation data, determine the fracture information in the vertical wellbore;

[0018] The fracture network model is obtained by fusing the main fracture trend surface and occurrence information of the horizontal section of the fractured well and the fracture information of the vertical wellbore into the initial fracture network model of the work area.

[0019] Preferably, a geomechanical model is constructed based on the fracture network model, well logging data, and drilling data, specifically including:

[0020] Based on the well-seismic fusion technology, the distribution trend of the structural surface in the work area is adjusted to obtain the structural model of the work area;

[0021] The vertical distribution characteristics of rock mechanics in a single well are characterized by well logging data. Based on the data obtained from laboratory tests and the vertical distribution characteristics of rock mechanics in a single well, dynamic and static transformations are completed through empirical formulas to obtain the work area attribute model.

[0022] The stress conditions of the model boundary were determined by trial and error, and the fracture network model, the structural model of the work area, and the property model of the work area were superimposed to obtain the geomechanical model of the work area.

[0023] Preferably, the direction cosine matrix is ​​calculated based on the angle between the normal vector of each fracture surface and the direction of the triaxial stress, specifically including:

[0024] The angle between the fracture surface normal and the direction of the maximum horizontal principal stress is calculated based on the angles between the fracture surface normal and the direction of the vertical principal stress. The specific formula is as follows:

[0025] ;

[0026] In the formula, The angle between the normal to the fracture surface and the direction of the minimum horizontal principal stress is given. The angle between the normal to the fracture surface and the direction of the maximum horizontal principal stress is given. The angle between the normal to the fracture surface and the direction of the vertical principal stress;

[0027] Based on the angles between the fracture surface normal and the direction of the minimum horizontal principal stress, and the angles between the fracture surface normal and the direction of the vertical principal stress, a direction cosine matrix is ​​constructed, with the following formula:

[0028] .

[0029] Preferably, the formula for calculating the activation risk of the fracture surface is:

[0030] ;

[0031] in,

[0032] ;

[0033] In the formula, The risk of activation of the fracture surface, Direction cosine matrix elements, , p The pore pressure inside the fracture. , and These are the maximum horizontal principal stress, the vertical principal stress, and the minimum horizontal principal stress, respectively.

[0034] Preferably, a system of linear algebraic equations is established with the shear displacement invariants of the three-dimensional triangular surface mesh elements as unknowns, and the three-dimensional slip distribution of the high-risk fracture surface is obtained by solving the equations. Specifically, this includes:

[0035] Multiple shear displacement discontinuities are placed on the high-risk fracture surface;

[0036] Based on the displacement discontinuity elements and the shear and normal boundary stresses, a system of linear equations for the continuous invariant shear displacement of each three-dimensional triangular surface mesh element is constructed, as follows:

[0037] ;

[0038] In the formula, , and It is the first Shear and normal boundary stresses on each element, with subscripts 1, 2, 3 indicating the three directions of the local coordinates; , and It is the first Discontinuous components of shear and normal displacement of each element; ( , =1,2,3) are the boundary influence coefficients.

[0039] Preferably, a linear hardening elastoplastic constitutive equation is introduced, using a three-dimensional slip distribution as the displacement boundary condition, and an eight-node hexahedral element is used to construct the casing mesh model. The nonlinear deformation of the casing is obtained through iterative calculation of the casing mesh model, and a casing inner diameter deformation contour map is plotted. Specifically, this includes:

[0040] The initial sleeve size model is divided into a sleeve mesh model using eight-node hexahedral elements;

[0041] The fracture shear slip is used as the displacement boundary condition of the casing model, and a linear hardening model is introduced to characterize the elastic-plastic properties of the casing.

[0042] The stress at the nodal points of the casing element is calculated using the stiffness matrix and Gaussian numerical integration. The stress and strain at the nodal points are also calculated. A three-dimensional cloud diagram of the casing deformation is plotted, and the deformation result curve is extracted to obtain the deformation of the inner diameter of the casing.

[0043] Preferably, the parameters of perforation cluster spacing, fracturing scale / displacement, and well spacing are dynamically optimized, specifically including:

[0044] After each iteration, the multi-point deformation cloud map of the casing is compared with the measured well logging data to calculate the error. If the error is greater than the set threshold, the perforation cluster spacing, fracturing scale, flow rate and well spacing are adjusted, and all actions of updating the fracture network model to the casing deformation calculation are re-executed. If the error is less than or equal to the set threshold, the current fracturing parameters are output as the optimal construction plan.

[0045] The above-mentioned at least one technical solution adopted in this invention can achieve the following beneficial effects:

[0046] To address the problem of overly idealistic calculations of casing deformation in shale reservoir fracturing techniques, this invention provides a method for calculating casing deformation in shale reservoir fracturing. This method utilizes a geomechanical model to construct a direction cosine matrix based on the magnitude and direction of the three principal stresses on the fracture surface. Activation risks are calculated, and high-risk fractures are screened out. The three-dimensional slip distribution is obtained through triangular mesh boundary element analysis. This establishes a multi-point casing deformation calculation method. The direction cosine matrix maps geometric angle information to the projection coefficients of the stress tensor, transforming the calculation of normal and shear stresses on the fracture surface from empirical estimation to tensor computation. Compared to the traditional simplification to vertical fractures, this method accurately reflects the slip potential of inclined fractures. Since the displacement discontinuity directly corresponds to fault slip, it avoids the accumulation of errors caused by secondary conversions of slip and deformation. This achieves dynamic feedback correction between the geological model and fracturing operation results, and completes dynamic optimization and adjustment of fracturing parameters, providing strong support for the efficient development of deep shale gas. Attached Figure Description

[0047] The accompanying drawings, which are included to provide a further understanding of this application and form part of this application, illustrate exemplary embodiments and are used to explain this application, but do not constitute an undue limitation of this application. In the drawings:

[0048] Figure 1 A flowchart illustrating a method for calculating the deformation of fracturing casing in shale reservoirs, provided by this invention;

[0049] Figure 2 A schematic diagram of the three-dimensional Mohr's circle stress state of a three-dimensional fault in spherical coordinates for a method of calculating the deformation of a fracturing casing in a shale reservoir provided by the present invention.

[0050] Figure 3 A schematic diagram of a three-dimensional displacement discontinuity model for a method of calculating the deformation of fracturing casing in shale reservoirs provided by the present invention;

[0051] Figure 4 The microseismic response and normalized fracture instability results of fracture-type reservoir fracturing in the A-platform well group provided by the present invention are based on a method for calculating casing deformation in shale reservoir fracturing.

[0052] Figure 5 The present invention provides a three-dimensional fracture fracturing instability shear slip result and slip curve of a method for calculating the deformation of fracturing casing in shale reservoirs.

[0053] Figure 6 The multi-arm caliber logging results and finite element method casing deformation results of well A1, which provide a method for calculating the deformation of fracturing casing in shale reservoirs according to the present invention;

[0054] Figure 7A schematic diagram of the fracture length under different perforation clusters in a single-stage fracturing operation, which is provided by the present invention for calculating the deformation of fracturing casing in shale reservoirs;

[0055] Figure 8 The fracturing results diagram of the fracture development, large approach angle fault and small approach angle fault development sections provided by the present invention are shown in the figure.

[0056] Figure 9 This invention provides a fracturing parameter optimization and control chart for a method of calculating the deformation of fracturing casing in shale reservoirs. Detailed Implementation

[0057] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this application will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments in the specification without creative effort are within the scope of protection of this application.

[0058] This embodiment presents a method for calculating casing deformation in shale reservoir fracturing. This method fully utilizes a large amount of microseismic data from the field to obtain the distribution characteristics of a fine fracture network in the reservoir, establishes multi-point casing deformation prediction in the wellbore, achieves dynamic feedback correction between the geological model and fracturing operation results, and completes dynamic optimization and adjustment of fracturing parameters, providing strong support for the efficient development of deep shale gas.

[0059] The technical solutions provided by the various embodiments of this application are described in detail below with reference to the accompanying drawings.

[0060] Figure 1 This is a schematic diagram of a method for calculating the deformation of fracturing casing in shale reservoirs according to the present invention, which specifically includes the following steps:

[0061] S101: Acquire microseismic data, imaging logging data, core observation data, and drilling and logging data generated during fracturing operations; construct a fracture network model based on the microseismic data, imaging logging data, and core observation data; construct a geomechanical model based on the fracture network model, logging data, and drilling and logging data; the geomechanical model includes triaxial stress data and directions for each fracture surface.

[0062] Optionally, a fracture network model is constructed based on microseismic data, imaging logging data, and core observation data. Specifically, this includes: determining the main fracture trend surface and attitude information of the horizontal section of the fractured well based on the microseismic data; the attitude information includes: fracture surface dip angle, fracture surface dip direction, and fracture surface area. An initial fracture network model for the work area is established based on the microseismic data, imaging logging data, and core observation data; the fracture information of the vertical wellbore is determined based on the imaging logging data and the core observation data; the main fracture trend surface and attitude information of the horizontal section of the fractured well and the fracture information of the vertical wellbore are integrated into the initial fracture network model for the work area to obtain the fracture network model.

[0063] Specifically, microseismic data from platform A's hydraulic fracturing were collected, and the b-values ​​of microseismic event points were calculated. The b-value is inversely proportional to stress concentration; higher b-values ​​observed in hydraulic fracturing are considered to represent a large number of natural fractures that open during high-pressure injection. Therefore, the b-value can be used to distinguish between matrix fracturing-induced event points and natural fracture propagation event points. The formula for calculating the b-value is as follows:

[0064] ;

[0065] Where: M is the magnitude; N is the number of earthquakes with a magnitude greater than or equal to M; a and b are constants reflecting seismic activity and seismotectonics. Event points representing natural fracture propagation are extracted from microseismic data. These event points are normalized and centered, and then, based on principal component analysis, they are clustered and fitted into multiple fracture surfaces to obtain the attitude information of the main fracture surfaces, specifically including dip angle, dip direction, and cross-sectional area. A fracture model is established based on seismic, imaging logging, and core observation results for the W work area. Large-scale fractures in the work area can be obtained through seismic data, while vertical wellbore fracture information can be obtained through imaging logging and core observation results. The natural fracture attitude information extracted from the microseismic data of platform A is added to the fracture model of the work area, updating and iterating the existing fracture model for the W work area to obtain a refined fracture model of platform A that accurately reflects the horizontal section of the fractured well.

[0066] S102: Based on the fracture network model, well logging data, and drilling data, construct a geomechanical model; the geomechanical model includes the magnitude and direction of the three principal stresses on each fracture surface.

[0067] Optionally, based on the fracture network model, well logging data, and drilling data, a geomechanical model is constructed, specifically including: adjusting the distribution trend of structural surfaces in the work area according to well-seismic fusion technology and geological understanding of the work area to obtain a structural model of the work area; characterizing the vertical distribution characteristics of rock mechanics in a single well through well logging data; and completing dynamic and static transformations through empirical formulas based on data obtained from laboratory tests and the vertical distribution characteristics of rock mechanics in a single well to obtain a property model of the work area; determining the stress conditions of the model boundary through trial and error, and superimposing the fracture network model, the structural model of the work area, and the property model of the work area to obtain a geomechanical model of the work area.

[0068] Specifically, seismic, well logging, and drilling data for work area W were collected. Well-seismic fusion technology was used to rationally adjust the structural surface distribution trend of work area and establish a structural model for work area W. Then, based on well logging data, the vertical distribution characteristics of rock mechanics in single wells were interpreted, and based on indoor experimental data, empirical formulas were used to complete the dynamic-static transformation, thereby constructing a property model for work area W. Finally, based on the trial-and-error method, the boundary stress conditions of the model were determined, a geomechanical model of work area W was established, and three-dimensional stress data of platform A, including the magnitude and direction of the three principal stresses, were obtained.

[0069] Specifically, the construction of the structural model for the work area includes: establishing a synthetic record based on well point sonic-density curves using seismic, well logging, and drilling data; establishing a precise depth-time correspondence through VSP calibration and strata control velocity field to achieve high-precision time-depth conversion; then tracking strata and identifying faults on the seismic profile, correcting the attitude using well fault points, and establishing a corner point grid framework; finally, using well points as hard constraints and seismic attributes as trend control, iteratively correcting the velocity field and fault morphology through geostatistical interpolation and Bayesian lithofacies fusion to achieve seamless unification of structural surface distribution trends and well-seismic data.

[0070] S103: Calculate the direction cosine matrix based on the angle between the normal of each fracture surface and the direction of the three principal stresses; calculate the activation risk value of each fracture surface based on the magnitude of the three principal stresses and the elements in the direction cosine matrix, and designate the fracture surfaces with activation risk values ​​greater than a preset threshold as high-risk fracture surfaces.

[0071] Optionally, based on the triaxial stress data and direction at the center point of each fracture surface, and combined with the dip angle of the fracture surface, the direction cosine matrix is ​​calculated, specifically including:

[0072] The angle between the fracture surface normal and the direction of the maximum horizontal principal stress is calculated based on the angles between the fracture surface normal and the direction of the vertical principal stress. The specific formula is as follows:

[0073] ;

[0074] In the formula, The angle between the normal to the fracture surface and the direction of the minimum horizontal principal stress is given. The angle between the normal to the fracture surface and the direction of the maximum horizontal principal stress is given. The angle between the normal to the fracture surface and the direction of the vertical principal stress;

[0075] Based on the angles between the fracture surface normal and the direction of the minimum horizontal principal stress, and the angles between the fracture surface normal and the direction of the vertical principal stress, a direction cosine matrix is ​​constructed, with the following formula:

[0076] .

[0077] Specifically, based on the geomechanical model of the work area, the triaxial stress characteristics at the fracture location are obtained. Combined with the fracture attitude and stress characteristics, the activation risk at different parts of the fracture surface can be determined. The distribution ranges from 0.6 to 1.0, and the fracture is in a state of mechanical activity. When the value is less than 0.6, the fracture is in a closed state.

[0078] Optionally, the formula for calculating the activation risk of the fracture surface is:

[0079] ;

[0080] in,

[0081] ;

[0082] In the formula, The risk of activation of the fracture surface, Direction cosine matrix elements, , p The pore pressure inside the fracture. , and These are the maximum horizontal principal stress, the vertical principal stress, and the minimum horizontal principal stress, respectively.

[0083] Specifically, based on the geomechanical results of platform A, combined with the microseismic data in step S101 ( Figure 4 The detailed fracture model of platform A, obtained from [previous data], was used to calculate the activation risk of each three-dimensional fracture on platform A using a fracture activation risk calculation model. Here, the maximum horizontal stress at platform A is 108 MPa, oriented 80° east of north; the minimum horizontal principal stress is 93.5 MPa, the vertical principal stress is 101.5 MPa, and the formation pressure is 80 MPa. Platform A experiences one major fracture (dip angle 90°, strike 39.5° east of north), with a calculated activation pressure increment of 7.72 MPa. A normalized fracture instability chart was also obtained; see [reference needed]. Figure 4 .

[0084] S104: The high-risk fracture surface is meshed into three-dimensional triangular surface mesh elements. A system of linear algebraic equations is established with the continuous invariant shear displacement of the three-dimensional triangular surface mesh elements as unknowns. The three-dimensional slip distribution of the high-risk fracture surface is obtained by solving the equations.

[0085] Optionally, a system of linear algebraic equations is established with the continuous invariant shear displacement of the three-dimensional triangular surface mesh element as unknowns, and the three-dimensional slip distribution of the high-risk fracture surface is obtained by solving the system. Specifically, this includes: placing multiple discontinuous shear displacement elements on the high-risk fracture surface; and constructing a system of linear equations for the continuous invariant shear displacement of each three-dimensional triangular surface mesh element based on the discontinuous displacement elements and the shear and normal boundary stresses, as follows:

[0086] ;

[0087] In the formula, , and It is the first Shear and normal boundary stresses on each element, with subscripts 1, 2, 3 indicating the three directions of the local coordinates; , and It is the first Discontinuous components of shear and normal displacement of each element; ( , =1,2,3) are the boundary influence coefficients.

[0088] Specifically, the fracture slip is solved using the three-dimensional displacement discontinuity method, which involves solving for the discontinuous elements in the calculation formula. The core elements include two components: stress and displacement. The final three-dimensional fracture slip is obtained by summing the various shear displacement elements.

[0089] Specifically, the displacement discontinuity method is used to calculate the three-dimensional fracture slip. In this method, the boundary is meshed into a three-dimensional triangular surface, where each face of the mesh acts as a triangular dislocation. The force diagram is shown below. Figure 2 The displacement of each fracture element is defined in the local coordinate system; see [link to relevant documentation]. Figure 3 The z-axis is along the normal displacement direction:

[0090] ;

[0091] To numerically implement a displacement discontinuity element, a displacement discontinuity D on the element is required. i The analytical solution. The general form of the displacement discontinuous element can be expressed as follows:

[0092] ;

[0093] Where G is the stiffness modulus; ν is Poisson's ratio; It's a kernel function.

[0094] ;

[0095] f x f xy f xyz Waiting is The partial derivatives are relative to the x, y, and z directions. Based on the solutions of the aforementioned constant three-dimensional displacement discontinuity element, a program for numerically solving representative boundary element problems can be developed. N displacement discontinuities are placed on the planar boundary. By considering the boundary conditions, a system of 3N linear algebraic equations concerning the unknown displacement discontinuities can be established.

[0096] Specifically, according to the body method, the problem of the original crack in a finite body can be divided into two subproblems: one is the external body without cracks, and the other is the crack. In the second subproblem, only the loads acting on the two surfaces of the crack are known. Then, all discontinuous displacement components can be solved by equation (2). Once the DDM algorithm model of three-dimensional fracture stress deformation is built, it can be directly solved and calculated. The specific steps are: create the fracture plane, define the boundary, the locking element of the boundary element and the interface; define the full space and elastic constants of the material, and create the data structure; define the boundary conditions, such as the stress, traction force, friction, etc. at the center of the element, and use the boundary element to calculate the slip distribution; draw the animation into a graph and extract the curve results.

[0097] Specifically, the high-risk fracture F on platform A was selected to calculate the fracture slip. The effective fracture surface length × height of fracture F is 226m × 20m, the boundary conditions are the aforementioned in-situ stress and formation pressure, and the material parameter is shale rock. The cloud map of the shear slip result of fracture F calculated based on the three-dimensional DDM method is shown below. Figure 5 The extraction result showed that the slip at the fracture center was 59.239 mm.

[0098] S105: Using the three-dimensional slip distribution as the displacement boundary condition, a linear hardening elastoplastic constitutive equation is introduced, and an eight-node hexahedral element is used to construct a casing mesh model; by iteratively calculating the casing mesh model, the nonlinear deformation of the casing is obtained, and a casing inner diameter deformation cloud map is drawn.

[0099] Optionally, a linear hardening elastoplastic constitutive equation is introduced, using the three-dimensional slip distribution as the displacement boundary condition, and an eight-node hexahedral element is used to construct the casing mesh model. The nonlinear deformation of the casing is obtained by iterative calculation of the casing mesh model, and the deformation cloud map of the casing inner diameter is plotted. Specifically, this includes: dividing the initial casing size model into a casing mesh model using eight-node hexahedral elements; using the fracture shear slip as the displacement boundary condition of the casing model, introducing a linear hardening model to characterize the elastoplastic characteristics of the casing; calculating the nodal stress of the casing element using the stiffness matrix and Gaussian numerical integration, calculating the nodal stress and strain, plotting the three-dimensional cloud map of the casing deformation, and extracting the deformation result curve to obtain the deformation of the casing inner diameter.

[0100] Specifically, the overall program design process is as follows: Given a casing size model, the casing mesh model is generated using C3D8 eight-node hexahedral elements; based on the fracture shear slip as the displacement boundary condition of the casing model, a linear hardening model is introduced to characterize the elasto-plastic characteristics of the casing; the Newton-Raphson iterative algorithm is used to solve the nonlinear deformation problem of the casing; the stiffness matrix and Gaussian numerical integral are used to calculate the stress at the casing element nodes; post-processing calculates the stress and strain at the nodes; the Patch method is used to realize the three-dimensional cloud map of the casing deformation; the deformation result curve is extracted to obtain the deformation of the casing inner diameter, specifically including:

[0101] (1) C3D8 element isoparametric transformation: To transform a geometrically regular element in the local coordinate system into a geometrically irregular element in the global coordinate system, a coordinate transformation needs to be established.

[0102] ;

[0103] The mechanical interpolation function for coordinates is in the form of:

[0104] ;

[0105] The interpolation function for displacement is in the form of:

[0106] ;

[0107] The mathematical process of substitution integration is as follows:

[0108] ;

[0109] (2) The formula for calculating the element stiffness matrix is ​​as follows:

[0110] ;

[0111] in, The strain-displacement matrix is ​​represented using global coordinates; The strain-displacement matrix is ​​expressed in isoparametric coordinates; The stress-strain matrix, It is a Jacobian matrix; , , , where is the Gaussian integral weight, and 2×2×2 Gaussian integral points are used.

[0112] (3) Three-dimensional unit Gaussian numerical integration:

[0113] ;

[0114] (4) Linear hardening elastoplastic constitutive equation:

[0115] ;

[0116] (5) The method of iteratively solving the nonlinear problem of the casing deformation process, which gradually approximates the minimum value of the objective function with the minimum value of the quadratic curve, has a fast convergence rate.

[0117] For a system of nonlinear equations, the matrix form is: The specific form is as follows:

[0118] ;

[0119] in, , Then the Jacobian matrix of the Frechet derivative of the function F(x) is expressed as:

[0120] ;

[0121] The iterative scheme for the nonlinear equation system is:

[0122] .

[0123] Specifically, based on the fracture slippage result obtained in step S102, considering the casing material parameters, cement sheath parameters, and the wellbore distance from the fracture center of 63.36m, the calculated casing slippage is 48.28mm. (See [link to relevant documentation]). Figure 6 The measured casing deformation during multi-arm wellbore logging was 52.72 mm, with an error of 9%. The casing used was TP140 (yield stress 965 MPa) steel grade. The material parameters for the wellbore model are shown in Table 1.

[0124] Table 1 Material parameters of the wellbore model

[0125]

[0126] In addition, the fracture network model and geomechanical model are iteratively updated based on newly added microseismic data. Combined with the prediction results of fracturing fracture length, stress disturbance range and casing deformation risk, the parameters of perforation cluster spacing, fracturing scale / discharge and well spacing are dynamically optimized.

[0127] Optionally, the parameters of perforation cluster spacing, fracturing scale / displacement, and well spacing can be dynamically optimized, specifically including:

[0128] After each iteration, the multi-point deformation cloud map of the casing is compared with the measured well logging data to calculate the error. If the error is greater than the set threshold, the perforation cluster spacing, fracturing scale, flow rate and well spacing are adjusted, and all actions of updating the fracture network model to the casing deformation calculation are re-executed. If the error is less than or equal to the set threshold, the current fracturing parameters are output as the optimal construction plan.

[0129] Specifically, as pressure operations continue, new microseismic data are constantly generated, allowing for multiple iterations to optimize the geomechanical model and the fine fracture model, and to complete the iterative optimization and upgrading of the risk model. Based on the length of the fracturing fracture and the range of stress disturbance, while taking into account the reservoir stimulation effect and the prevention of fracturing / cascade transformation accidents, dynamic optimization and adjustment of fracturing parameters can be achieved, mainly including (perforation cluster spacing optimization, fracturing scale / displacement optimization, and well spacing). Figure 7 It is the fracture length under single-stage fracturing conditions of 4 and 8 clusters. Figure 8 It is the result of the propagation of hydraulic fracturing fractures in naturally fractured and fracture-developed sections. Figure 9 The chart for optimizing and controlling fracturing parameters of Platform A is designed based on the detailed fracture model and stress model of Platform A, and the optimization design of fracturing parameters of the platform is carried out, while preventing the occurrence of casing deformation.

[0130] This patented technology establishes a fully three-dimensional realistic fracture model (with dip and strike) based on microseismic fracture characterization technology. Based on this, a calculation method integrating a three-dimensional activation risk model, a three-dimensional fracture slip model, and a three-dimensional casing stress-deformation model under realistic three-dimensional fracture conditions is constructed, enabling rapid solution of realistic casing variables. Furthermore, it organically combines fracturing operations with casing deformation risk prediction, completing iterative optimization and upgrading of the risk model, achieving dynamic optimization and adjustment of fracturing parameters, and enabling subsequent fracturing operation decisions.

[0131] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this invention.

Claims

1. A method for calculating the deformation of fracturing casing in shale reservoirs, characterized in that, include: Acquire microseismic data, imaging logging data, core observation data, and drilling data generated in shale reservoirs during fracturing operations; A fracture network model was constructed based on microseismic data, imaging logging data, and core observation data. Based on the fracture network model, well logging data, and drilling data, a geomechanical model is constructed; the geomechanical model includes the magnitude and direction of the three principal stresses on each fracture surface; Calculate the direction cosine matrix based on the angle between the normal of each fracture surface and the direction of the three principal stresses; calculate the activation risk value of each fracture surface based on the magnitude of the three principal stresses and the elements in the direction cosine matrix, and designate the fracture surfaces with activation risk values ​​greater than a preset threshold as high-risk fracture surfaces. The high-risk fracture surface is meshed into three-dimensional triangular surface mesh elements. A system of linear algebraic equations is established with the continuous invariant shear displacement of the three-dimensional triangular surface mesh elements as unknowns. The three-dimensional slip distribution of the high-risk fracture surface is obtained by solving the system. Specifically, this includes: placing multiple shear displacement discontinuities on the high-risk fracture surface; and constructing a system of linear equations for the continuous invariant shear displacement of each three-dimensional triangular surface mesh element based on the displacement discontinuities and the shear and normal boundary stresses. The formula is as follows: ； In the formula, , and It is the first i Shear and normal boundary stresses on each element, with subscripts 1, 2, 3 indicating the three directions of the local coordinates; , and It is the first j Discontinuous components of shear and normal displacement of each element; ( m , n =1,2,3) are the boundary influence coefficients; Using the three-dimensional slip distribution as the displacement boundary condition, a linear hardening elastoplastic constitutive equation is introduced, and an eight-node hexahedral element is used to construct a casing mesh model. The nonlinear deformation of the casing is obtained through iterative calculation of the casing mesh model, and a casing inner diameter deformation cloud map is plotted. Specifically, this includes: dividing the initial casing size model into a casing mesh model using eight-node hexahedral elements; using the fracture shear slip as the displacement boundary condition of the casing model, introducing a linear hardening model to characterize the elastoplastic properties of the casing; calculating the nodal stress of the casing elements using the stiffness matrix and Gaussian numerical integration, calculating the nodal stress-strain, plotting the three-dimensional casing deformation cloud map, and extracting the deformation result curve to obtain the casing inner diameter deformation.

2. The method for calculating the deformation of fracturing casing in shale reservoirs as described in claim 1, characterized in that, The method also includes: iteratively updating the fracture network model and geomechanical model based on newly added microseismic data, and dynamically optimizing the parameters of perforation cluster spacing, fracturing scale / discharge and well spacing by combining the fracturing fracture length, stress disturbance range and casing deformation risk prediction results.

3. The method for calculating the deformation of fracturing casing in shale reservoirs as described in claim 1, characterized in that, The construction of a fracture network model based on microseismic data, imaging logging data, and core observation data specifically includes: Based on microseismic data, the main fracture trend surface and attitude information of the horizontal section of the fractured well are determined; the attitude information includes: fracture surface dip angle, fracture surface dip direction, and fracture surface area. An initial fracture network model for the work area was established based on microseismic data, imaging logging data, and core observation data. Based on the imaging logging data and the core observation data, determine the fracture information in the vertical wellbore; The fracture network model is obtained by fusing the main fracture trend surface and occurrence information of the horizontal section of the fractured well and the fracture information of the vertical wellbore into the initial fracture network model of the work area.

4. The method for calculating the deformation of fracturing casing in shale reservoirs as described in claim 1, characterized in that, The construction of the geomechanical model based on the fracture network model, well logging data, and drilling data specifically includes: Based on the well-seismic fusion technology, the distribution trend of the structural surface in the work area is adjusted to obtain the structural model of the work area; The vertical distribution characteristics of rock mechanics in a single well are characterized by well logging data. Based on the data obtained from laboratory tests and the vertical distribution characteristics of rock mechanics in a single well, dynamic and static transformations are completed through empirical formulas to obtain the work area attribute model. The stress conditions of the model boundary were determined by trial and error, and the fracture network model, the structural model of the work area, and the property model of the work area were superimposed to obtain the geomechanical model of the work area.

5. The method for calculating the deformation of fracturing casing in shale reservoirs as described in claim 1, characterized in that, The calculation of the direction cosine matrix based on the angle between the normal vector of each fracture surface and the direction of the triaxial stress specifically includes: The angle between the fracture surface normal and the direction of the maximum horizontal principal stress is calculated based on the angles between the fracture surface normal and the direction of the vertical principal stress. The specific formula is as follows: ； In the formula, The angle between the normal to the fracture surface and the direction of the minimum horizontal principal stress is given. The angle between the normal to the fracture surface and the direction of the maximum horizontal principal stress is given. The angle between the normal to the fracture surface and the direction of the vertical principal stress; Based on the angles between the fracture surface normal and the direction of the minimum horizontal principal stress, and the angles between the fracture surface normal and the direction of the vertical principal stress, a direction cosine matrix is ​​constructed, with the following formula: 。 6. The method for calculating the deformation of fracturing casing in shale reservoirs as described in claim 5, characterized in that, The formula for calculating the activation risk of the fracture surface is: ； in, ； In the formula, The risk of activation of the fracture surface, Direction cosine matrix elements, , p The pore pressure inside the fracture. , and These are the maximum horizontal principal stress, the vertical principal stress, and the minimum horizontal principal stress, respectively.

7. The method for calculating the deformation of fracturing casing in shale reservoirs as described in claim 2, characterized in that, The dynamically optimized parameters for perforation cluster spacing, fracturing scale / displacement, and well spacing specifically include: After each iteration, the multi-point deformation cloud map of the casing is compared with the measured well logging data to calculate the error. If the error is greater than the set threshold, the perforation cluster spacing, fracturing scale, flow rate and well spacing are adjusted, and all actions of updating the fracture network model to the casing deformation calculation are re-executed. If the error is less than or equal to the set threshold, the current fracturing parameters are output as the optimal construction plan.

Images