Shale reservoir fracturing casing deformation calculation method

Abstract

The invention discloses a shale reservoir fracturing casing deformation calculation method, and relates to the technical field of shale oil and gas development. The method comprises the following steps: fusing microseism, logging, rock core and drilling and logging data, and constructing a fracture network-geomechanics coupling model; calculating a direction cosine matrix according to the fracture three-dimensional stress and the inclination angle, evaluating an activation risk, and screening out high-risk fractures; the triangular mesh boundary element solves the slippage as a displacement boundary, and an eight-node hexahedron elastic-plastic casing mesh is utilized to iteratively calculate a casing variable and draw a cloud picture; and a micro-seismic data real-time iteration model is newly added along with construction, and the perforation cluster distance, the fracturing scale, the displacement and the well spacing are dynamically optimized. According to the method, a large amount of on-site microseismic data is fully utilized, reservoir fracture fine description and slippage risk analysis are achieved, shaft multi-point casing variable prediction is completed, and key technical support is provided for deep shale gas development fracturing scheme decision and risk avoidance.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of shale oil and gas development, and particularly relates to a shale reservoir fracturing casing deformation amount calculation method. BACKGROUND

[0002] In the fracturing operation of deep shale gas, a large amount of casing deformation (hereinafter referred to as casing deformation) and even channeling phenomenon are caused by complex geological conditions and fracturing parameters. The casing deformation problem affects the wellbore integrity, causes the combined pressure loss of the horizontal well section to occur frequently, and even some wells directly lose productivity. The channeling phenomenon is more difficult, which causes a large amount of adjacent well productivity to be reduced, seriously affects the fracturing construction process, and affects the production rhythm. Existing research shows that the main factor causing casing deformation is that fluid causes fault shear slip, and the fault shear acts on the casing to cause casing deformation. According to the deformation degree of the casing inner diameter, the casing deformation is divided into five levels in the field, and corresponding fracturing reconstruction schemes can be adopted according to the casing deformation level to reduce the influence of casing deformation. Therefore, accurate calculation and prediction of the casing deformation amount are crucial for deep shale gas development.

[0003] In the prior art, the fracture morphology is considered by using idealized two-dimensional vertical fractures, and the two-dimensional Mohr circle is used to solve the activation risk and slip amount. However, the prediction result of the casing deformation amount of the shale reservoir in the prior art is too idealized and not accurate enough, so that it cannot truly reflect the casing deformation risk.

[0004] Therefore, there is an urgent need for a shale reservoir fracturing casing deformation amount calculation method to effectively improve the prediction accuracy of the casing deformation amount of the rock reservoir. SUMMARY

[0005] Therefore, it is necessary to provide a shale reservoir fracturing casing deformation amount calculation method in view of the above technical problems.

[0006] The present application adopts the following technical scheme: The present application provides a shale reservoir fracturing casing deformation amount calculation method, comprising: Obtaining microseismic data, imaging logging data, core observation data and drilling and logging data generated by the shale reservoir under the condition of fracturing construction; constructing a fracture network model according to the microseismic data, the imaging logging data and the core observation data; Based on the fracture network model, the logging data and the drilling and logging data, a geomechanical model is constructed; the geomechanical model includes the size and direction of the three principal stresses of each fracture surface; According to the included angle between the normal line of each fracture surface and the direction of the three principal stresses, a direction cosine matrix is calculated; according to the size of the three principal stresses and the elements in the direction cosine matrix, the activation risk value of each fracture surface is calculated, and the fracture surface corresponding to the activation risk value greater than the preset threshold value is taken as a high-risk fracture surface; The high-risk fracture surface is meshed into three-dimensional triangular curved surface grid units, linear algebraic equations are established by taking shear displacement continuity invariants of the three-dimensional triangular curved surface grid units as unknowns, and three-dimensional slip distribution of the high-risk fracture surface is obtained by solving; A linear hardening elastic-plastic constitutive equation is introduced by taking the three-dimensional slip distribution as a displacement boundary condition, and a casing grid model is constructed by using eight-node hexahedral elements; and a casing nonlinear deformation is obtained by iteratively calculating the casing grid model, and a casing inner diameter deformation cloud chart is drawn.

[0007] Preferably, the method further comprises: iteratively updating the fracture network model and the geomechanical model according to the new microseismic data, dynamically optimizing the perforation cluster spacing, the fracturing scale / displacement and the well spacing parameters in combination with the fracturing fracture length, the stress disturbance range and the casing deformation risk prediction results.

[0008] Preferably, the fracture network model is constructed according to the microseismic data, the imaging logging data and the core observation data, and specifically comprises: The fracture network model is constructed according to the microseismic data, the imaging logging data and the core observation data, and specifically comprises: The fracture network model is constructed according to the microseismic data, the imaging logging data and the core observation data, and specifically comprises: The fracture network model is constructed according to the microseismic data, the imaging logging data and the core observation data, and specifically comprises: The fracture network model is constructed according to the microseismic data, the imaging logging data and the core observation data, and specifically comprises:

[0009] Preferably, the geomechanical model is constructed based on the fracture network model, the logging data and the drilling and logging data, and specifically comprises: The structure model of the work area is obtained by adjusting the structure surface distribution trend of the work area according to the well-to-seismic fusion technology; The attribute model of the work area is obtained by completing dynamic and static conversion through an empirical formula according to the data obtained from the laboratory test and the longitudinal distribution characteristics of the single well rock mechanics, and the attribute model of the work area is obtained by completing dynamic and static conversion through an empirical formula according to the data obtained from the laboratory test and the longitudinal distribution characteristics of the single well rock mechanics; The stress condition of the model boundary is determined by trial and error, and the fracture network model, the structure model of the work area and the attribute model of the work area are superimposed to obtain the geomechanical model of the work area.

[0010] Preferably, the direction cosine matrix is calculated according to the angle between the normal vector of each fracture surface and the direction of the three-dimensional stress, and specifically comprises: The angle between the normal line of the fracture surface and the direction of the minimum horizontal principal stress is calculated according to the angle between the normal line of the fracture surface and the direction of the maximum horizontal principal stress and the angle between the normal line of the fracture surface and the direction of the vertical principal stress, and the specific formula is: ; wherein, is the angle between the normal of the fracture surface and the direction of the minimum horizontal principal stress, is the angle between the normal of the fracture surface and the direction of the maximum horizontal principal stress, is the angle between the normal of the fracture surface and the direction of the vertical principal stress; According to the angle between the normal of the fracture surface and the direction of the minimum horizontal principal stress and the angle between the normal of the fracture surface and the direction of the vertical principal stress, a direction cosine matrix is constructed, and the formula is: .

[0011] Preferably, the formula for calculating the activation risk of the fracture surface is: ; wherein, ; wherein, is the activation risk of the fracture surface, is the element of the direction cosine matrix , , p is the pore pressure inside the fracture, , and are the maximum horizontal principal stress, the vertical principal stress and the minimum horizontal principal stress, respectively.

[0012] Preferably, a linear algebraic equation group is established with the shear displacement continuous invariant of the three-dimensional triangular curved surface grid element as the unknown, and the three-dimensional slip distribution of the high-risk fracture surface is solved, specifically including: a plurality of shear displacement discontinuity elements are placed on the high-risk fracture surface; According to the displacement discontinuity element and the shear and normal boundary stress, a linear equation group of the shear displacement continuous invariant of each three-dimensional triangular curved surface grid element is constructed, and the formula is: ; wherein, , and are the shear and normal boundary stresses on the first element, and the subscripts 1, 2 and 3 represent three directions of local coordinates; , and are the shear and normal displacement discontinuity components of the first element; , =1, 2, 3) are boundary influence coefficients.

[0013] ​Preferably, the linear hardening elastic-plastic constitutive equation is introduced with the three-dimensional slip amount distribution as the displacement boundary condition, and an eight-node hexahedral element is adopted to construct the casing grid model; through iterative calculation on the casing grid model, the nonlinear deformation amount of the casing is obtained, and a casing inner diameter deformation cloud chart is drawn, specifically including: The initial casing size model is divided into a casing grid model through an eight-node hexahedral element; The fracture shear slip amount is taken as the displacement boundary condition of the casing model, and a linear hardening model is introduced to represent the elastic-plastic characteristics of the casing. The stress of the casing element node is calculated by using the stiffness matrix and the Gauss numerical integral, and the node stress strain is calculated, a three-dimensional cloud chart of the casing is drawn, and the casing inner diameter deformation amount is obtained by extracting the deformation amount result curve.

[0014] Preferably, the perforation cluster spacing, the fracturing scale / displacement and the well spacing parameters are dynamically optimized, specifically including: After each iteration, the casing multi-point deformation cloud chart is compared with the measured caliper logging data, and the error is calculated; if the error is greater than a set threshold, the perforation cluster spacing, the fracturing scale, the displacement and the well spacing are adjusted, and the whole action of the fracture network model updating to the casing deformation amount calculation is re-executed; if the error is less than or equal to the set threshold, the current fracturing parameters are output as the optimal construction scheme.

[0015] The above at least one technical scheme adopted by the present application can achieve the following beneficial effects: In view of the technical problem that the calculation result of the shale reservoir fracturing casing deformation amount in the prior art is too idealistic, the present application provides a shale reservoir fracturing casing deformation amount calculation method, the size and direction of the three-direction principal stress of the fracture surface are obtained through the constructed geomechanical model, the direction cosine matrix is constructed and the activation risk is calculated and the high-risk fracture is screened out, the three-dimensional slip distribution is obtained through the triangular grid boundary element solution, the calculation of the casing deformation amount at the wellbore multi-point position is established, the geometric angle information is mapped into the projection coefficient of the stress tensor by the direction cosine matrix, the normal stress and the shear stress of the fracture surface are calculated from the empirical estimation to the tensor operation, compared with the traditional method of simplifying the vertical fracture, the slip potential of the inclined fracture can be truly reflected, since the displacement discontinuity directly corresponds to the fault slip, the error accumulation caused by the secondary conversion of the slip and deformation is avoided, the dynamic mutual feedback correction between the geologic model and the fracturing construction result is realized, the dynamic optimization adjustment of the fracturing parameters is completed, and strong support is provided for efficient development of deep shale gas. BRIEF DESCRIPTION OF DRAWINGS

[0016] The drawings described herein are used to provide further understanding of the present application, and form a part of the present application, the illustrative embodiments of the present application and the description thereof are used to explain the present application, and do not constitute improper limitations on the present application. In the drawings: Figure 1A flowchart of a shale reservoir fracturing casing deformation calculation method provided by the present application is provided. Figure 2 A three-dimensional fault in the spherical coordinate system under the three-dimensional Mohr circle stress state diagram of the shale reservoir fracturing casing deformation calculation method provided by the present application is provided. Figure 3 A three-dimensional displacement discontinuity model diagram of the shale reservoir fracturing casing deformation calculation method provided by the present application is provided. Figure 4 A platform well group fracturing microseismic response and normalized fracture instability result of the shale reservoir fracturing casing deformation calculation method provided by the present application is provided. Figure 5 A three-dimensional fracture fracturing instability shear slip result and slip amount curve diagram of the shale reservoir fracturing casing deformation calculation method provided by the present application is provided. Figure 6 A platform 1 well multi-arm caliper logging result and finite element method casing deformation result diagram of the shale reservoir fracturing casing deformation calculation method provided by the present application is provided. Figure 7 A single-stage fracturing different perforation cluster under fracturing fracture length diagram of the shale reservoir fracturing casing deformation calculation method provided by the present application is provided. Figure 8 A fracture development, large approach angle fault and small approach angle fault development layer fracturing result diagram of the shale reservoir fracturing casing deformation calculation method provided by the present application is provided. Figure 9 A fracturing parameter optimization control chart of the shale reservoir fracturing casing deformation calculation method provided by the present application is provided. DETAILED DESCRIPTION

[0017] In order to make the purpose, technical scheme and advantages of the present application clearer, the technical scheme of the present application will be described clearly and completely below in combination with specific embodiments of the present application and corresponding drawings. Obviously, the described embodiments are only some of the embodiments of the present application, not all the embodiments. Based on the embodiments in the specification, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application.

[0018] The present application provides a shale reservoir fracturing casing deformation calculation method. This method makes full use of a large amount of microseismic data on site, obtains fine fracture network distribution characteristics of the reservoir, establishes multi-point casing deformation variable prediction of the wellbore, realizes dynamic feedback correction of the geological model and fracturing construction results, completes dynamic optimization and adjustment of the fracturing parameters, and provides strong support for efficient development of deep shale gas. The technical scheme provided by each embodiment of the present application will be described in detail below in combination with the drawings.

[0019] Figure 1 It is a shale reservoir fracturing casing deformation calculation method flowchart, specifically comprising the following steps: S101: Obtain microseismic data, imaging logging data, core observation data and drilling and logging data generated by fracturing operation; construct a fracture network model according to 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 three-dimensional stress data and direction of each fracture surface.

[0020] Optionally, the fracture network model is constructed according to the microseismic data, imaging logging data and core observation data, specifically comprising: determining the main fracture trend surface and occurrence information of the fracture in the horizontal section of the fracturing well according to the microseismic data; the occurrence information includes: fracture surface dip angle, fracture surface trend and fracture surface area. Establishing an initial fracture network model of the work area according to the microseismic data, imaging logging data and core observation data; determining the straight well borehole fracture information according to the imaging logging data and the core observation data; fusing the main fracture trend surface and occurrence information of the fracture in the horizontal section of the fracturing well and the straight well borehole fracture information to the initial fracture network model of the work area to obtain the fracture network model.

[0021] Specifically, collect the fracturing microseismic data of platform A, calculate the b value of the microseismic event point, the b value is inversely proportional to stress concentration, and the higher b value observed in hydraulic fracturing is considered to represent a large number of natural fractures opened during high-pressure injection. Therefore, the b value can be used to distinguish between matrix fracture-induced event points and natural fracture expansion event points. The b value calculation formula is as follows: ; Where: M is the magnitude; N is the number of earthquakes with magnitude greater than or equal to M; a and b are constants, which reflect seismic activity and seismic structure; extract natural fracture expansion event points from microseismic data, normalize and centralize the data of these event points, then based on principal component analysis, cluster and fit the event points into multiple fracture surfaces, and obtain the occurrence information of the main fracture surface, including dip angle, trend and fracture surface area; based on the seismic, imaging logging and core observation results of work area W, a fracture model is established, wherein the large-scale fractures in the work area can be obtained through seismic data, and the straight well borehole fracture information can be obtained through imaging logging and core observation results. Add the natural fracture occurrence information extracted by the microseismic of platform A to the fracture model of the work area, update and iterate the established fracture model of work area W, and obtain the accurate A platform fine fracture model that can accurately reflect the horizontal section of the fracturing well.

[0022] S102: Construct a geomechanical model based on the fracture network model, logging data and drilling and logging data; the geomechanical model includes the size and direction of the three principal stresses of each fracture surface.

[0023] Optionally, based on the fracture network model, logging data and drilling and logging data, a geomechanical model is constructed, specifically including: adjusting the distribution trend of the structural surface of the work area according to the well-seismic fusion technology and the geological understanding of the work area to obtain a structural model of the work area; representing the vertical distribution characteristics of single well rock mechanics through logging data, and completing dynamic and static conversion through empirical formula according to the data obtained from indoor test and the vertical distribution characteristics of single well rock mechanics to obtain an attribute model of the work area; determining the stress condition of the model boundary through trial and error, and superimposing the fracture network model, the structural model of the work area and the attribute model of the work area to obtain a geomechanical model of the work area.

[0024] Specifically, seismic, logging, drilling and logging data of the W work area are collected, well-seismic fusion technology is adopted, the distribution trend of the structural surface of the work area is reasonably adjusted, and a structural model of the W work area is established; then the vertical distribution characteristics of single well rock mechanics are explained according to logging data, and dynamic and static conversion is completed through empirical formula based on indoor test data, and then an attribute model of the W work area is constructed; finally, the stress condition of the model boundary is determined based on trial and error, a geomechanical model of the W work area is established, and three-dimensional stress data of the A platform are obtained, including the size and direction of the three principal stresses.

[0025] Specifically, the construction of the structural model of the work area includes: generating a composite record based on the well point acoustic wave-density curve established based on the seismic, logging, drilling and logging data of the work area, and establishing a depth-time accurate correspondence through VSP calibration and layer-controlled velocity field to complete high-precision time-depth conversion; then the horizon is tracked on the seismic profile, the faults are identified, and the occurrence is corrected with the well breakpoint to establish an angle point grid framework; finally, the well point is taken as a hard constraint, the seismic attribute is taken as a trend control, the velocity field and the fault shape are iteratively corrected through geostatistical interpolation and Bayes facies fusion, and the seamless integration of the structural surface distribution trend and the well-seismic data is realized.

[0026] S103: Calculate the direction cosine matrix according to the included angle between the normal line of each fracture surface and the direction of the three principal stresses; calculate the activation risk value of each fracture surface according to the size of the three principal stresses and the elements in the direction cosine matrix, and take the fracture surface corresponding to the activation risk value greater than the preset threshold value as the high-risk fracture surface.

[0027] Optionally, the direction cosine matrix is calculated according to the three-dimensional stress data and the direction of the center point of each fracture surface in combination with the dip angle of the fracture surface, specifically including: The included angle between the normal line of the fracture surface and the direction of the minimum horizontal principal stress is calculated according to the included angle between the normal line of the fracture surface and the direction of the maximum horizontal principal stress and the included angle between the normal line of the fracture surface and the direction of the vertical principal stress, and the specific formula is: ; wherein, is the angle between the normal of the fracture surface and the direction of the minimum horizontal principal stress, is the angle between the normal of the fracture surface and the direction of the maximum horizontal principal stress, is the angle between the normal of the fracture surface and the direction of the vertical principal stress; According to the angle between the normal of the fracture surface and the direction of the minimum horizontal principal stress and the angle between the normal of the fracture surface and the direction of the vertical principal stress, a direction cosine matrix is constructed, and the formula is: .

[0028] Specifically, based on the geomechanical model of the work area, the three-dimensional stress characteristics of the location of the fracture are obtained, and the activation risk of different parts of the fracture surface can be obtained in combination with the fracture occurrence and stress characteristics, when is distributed in 0.6-1.0, the fracture is in a mechanical activity state, when is less than 0.6, the fracture is in a closed state.

[0029] Optionally, the formula for calculating the activation risk of the fracture surface is: ; wherein, ; wherein, is the activation risk of the fracture surface, is an element of the direction cosine matrix , , p is the pore pressure inside the fracture, , and are the maximum horizontal principal stress, the vertical principal stress and the minimum horizontal principal stress, respectively.

[0030] Specifically, based on the geomechanical results of the A platform, in combination with the fine fracture model of the A platform obtained in step S101 (microseismic data Figure 4 ), the activation risk of each three-dimensional fracture of the A platform can be obtained by using the fracture activation risk calculation model. Here, the maximum horizontal stress of the A platform is 108 MPa, the direction is 80° north of east; the minimum horizontal principal stress is 93.5 MPa, the vertical principal stress is 101.5 MPa, and the formation pressure is 80 MPa. The A platform experiences one large fracture (dip angle 90°, strike 39.5° north of east), the calculated activation pressure increment is 7.72 MPa, and the normalized fracture instability chart is obtained, see Figure 4 .

[0031] S104: mesh the high-risk fracture surface into three-dimensional triangular curved surface grid units, establish a linear algebraic equation group with the shear displacement continuous invariant of the three-dimensional triangular curved surface grid unit as the unknown, and solve to obtain the three-dimensional slip distribution of the high-risk fracture surface.

[0032] Optionally, the three-dimensional slip distribution of the high-risk fracture surface is solved by establishing a linear algebraic equation group with the shear displacement continuous invariant of the three-dimensional triangular curved surface grid unit as the unknown, specifically including: placing a plurality of shear displacement discontinuity elements on the high-risk fracture surface; constructing a linear equation group of the shear displacement continuous invariant of each three-dimensional triangular curved surface grid unit according to the displacement discontinuity element and the shear and normal boundary stress, the formula being: ; In the formula, , and are the shear and normal boundary stresses on the first unit, and the subscripts 1, 2 and 3 represent three directions of local coordinates; , and are the shear and normal displacement discontinuity components of the first unit; , =1, 2, 3) are boundary influence coefficients.

[0033] Specifically, the fracture slip is solved by the three-dimensional displacement discontinuity method, that is, the discontinuity element in the formula is calculated, and the core elements include two, one is stress and the other is displacement. Summing up each shear displacement element can obtain the final three-dimensional fracture slip.

[0034] Specifically, the three-dimensional fracture slip is calculated by the displacement discontinuity method, in which the boundary is meshed into a three-dimensional triangular curved surface, and each face of the mesh acts as a triangular dislocation, and the force diagram is shown in Figure 2 . The displacement of each fracture unit is defined in the local coordinate system, as shown in Figure 3 , and the z axis is along the normal displacement direction: ; In order to realize the displacement discontinuity unit numerically, the analytical solution of the displacement discontinuity D i on the unit is needed. The general form of the displacement discontinuity unit can be expressed as follows, ; where G is the stiffness modulus; v is the Poisson's ratio; and is the kernel function, ; f​x , f xy , f xyz etc. are With respect to the partial derivatives in the x, y and z directions. Based on the solution of the constant three-dimensional displacement discontinuity element described above, a program for numerically solving the representative boundary element problem can be developed. N displacement discontinuity elements are placed on the plane boundary. By considering the boundary conditions, a 3N linear algebraic equation system about the unknown displacement discontinuity components can be established.

[0035] Specifically, according to the body force method, the original crack problem in the finite body can be divided into two sub-problems: one is the outer body without cracks, and the other is the crack. In the second sub-problem, the load acting only on the two surfaces of the crack is known. Then, all displacement discontinuity components can be solved by equation (2). When the three-dimensional DDM algorithm model of the fracture stress deformation is built, it can be directly calculated. The specific steps are: creating fault surface, defining boundary, locking element and interface of boundary element; defining the full space and elastic constant of the material, creating data structure; defining boundary conditions, stress, traction, friction, etc. at the center of the element, using boundary element to calculate the slip distribution; drawing the graph and extracting the curve result.

[0036] Specifically, the high-risk fracture F on platform A is selected to calculate the fracture slip. The effective fracture surface length x height of the F fracture is 226m x 20m, the boundary condition is the above-mentioned ground stress and formation pressure, and the material parameter is shale rock. The shear slip result cloud map of the F fracture calculated based on the three-dimensional DDM method is shown in Figure 5 , and the extraction result of the fracture center slip is 59.239mm.

[0037] S105: taking the three-dimensional slip distribution as the displacement boundary condition, introducing a linear hardening elastic-plastic constitutive equation, and constructing a casing grid model using eight-node hexahedral elements; by iterative calculation on the casing grid model, the nonlinear deformation of the casing is obtained, and the casing inner diameter deformation cloud map is drawn.

[0038] Alternatively, taking the three-dimensional slip distribution as the displacement boundary condition, introducing a linear hardening elastic-plastic constitutive equation, and constructing a casing grid model using eight-node hexahedral elements; by iterative calculation on the casing grid model, the nonlinear deformation of the casing is obtained, and the casing inner diameter deformation cloud map is drawn. Specifically, the initial casing size model is divided into a casing grid model by eight-node hexahedral elements; the fracture shear slip is taken as the displacement boundary condition of the casing model, and a linear hardening model is introduced to represent the elastic-plastic characteristics of the casing; the stress of the casing element node is calculated by using the stiffness matrix and Gauss numerical integration, and the node stress and strain are calculated, and the three-dimensional cloud map of the casing is drawn, and the deformation result curve is obtained to obtain the casing inner diameter deformation.

[0039] Specifically, the overall program design specific process is: given the casing size model, the casing grid model is divided by using C3D8 eight-node hexahedral element; the linear hardening model is introduced to represent the elastic-plastic characteristics of the casing according to the fracture shear slip amount as the displacement boundary condition of the casing model, the Newton-Raphson iteration algorithm is used to solve the nonlinear deformation problem of the casing, and the stiffness matrix and the Gaussian numerical integration are used to calculate the stress of the casing element node; the post-processing calculation node stress and strain are carried out, the Patch method is used to realize the drawing of the casing three-dimensional cloud picture, and the casing inner diameter deformation is obtained by extracting the deformation result curve, which specifically includes: (1) C3D8 element isoparametric transformation: the element with regular geometric shape in local coordinates is converted into the element with irregular geometric shape in the global coordinate system, and a coordinate transformation needs to be established.

[0040] ; The mechanical interpolation function form of the coordinate is: ; The interpolation function form of the displacement is: ; The mathematical form of the element stiffness matrix is as follows: ; (2) Calculation of element stiffness matrix, the formula is: ; Wherein, is the strain-displacement matrix expressed by the global coordinate; is the strain-displacement matrix expressed by the isoparametric coordinate; is the stress-strain matrix, is the Jacobian matrix; , , , and

[0041] (3) Three-dimensional element Gaussian numerical integration: ; (4) Linear hardening elastic-plastic constitutive equation: ; (5) Iterative solution is used to calculate the nonlinear problem of the casing deformation process, and the minimum value of the quadratic curve gradually approaches the minimum value of the objective function, so the convergence rate is fast.

[0042] For the nonlinear equation set, the matrix form is , and the specific form is: ; Wherein, , . Then the Jacobian matrix of the Frechet derivative of the function F(x) is given by: ; Then the iterative scheme for the nonlinear system is .

[0043] Specifically, based on the fracture slip amount result obtained in step S102, considering the casing material parameters, cement sheath parameters and the distance of 63.36 m from the fracture center to the wellbore, the casing slip amount is calculated to be 48.28 mm, see Figure 6 ; the casing variable measured by the multi-arm caliper logging is 52.72 mm, with an error of 9%. The casing used is TP140 (yield stress 965 MPa) steel grade casing, and the wellbore model material parameter values are shown in Table 1:

[0044] Table 1 Wellbore model material parameters

[0045] In addition, the fracture network model and the geomechanical model are iteratively updated according to the newly added microseismic data, and the perforation cluster spacing, fracturing scale / displacement and well spacing parameters are dynamically optimized in combination with the fracturing fracture length, stress disturbance range and casing deformation risk prediction results.

[0046] Optionally, the perforation cluster spacing, fracturing scale / displacement and well spacing parameters are dynamically optimized, specifically including: After each iteration, the casing multi-point deformation cloud map is compared with the measured caliper logging data, and the error is calculated; if the error is greater than a set threshold, the perforation cluster spacing, fracturing scale, displacement and well spacing are adjusted, and the entire action of fracture network model updating to casing deformation calculation is re-executed; if the error is less than or equal to the set threshold, the current fracturing parameters are output as the optimal construction scheme.

[0047] Specifically, since the pressure operation is continuously carried out, new microseismic data can be continuously generated, and the geomechanical model and the fine fracture model can be iteratively optimized to complete the iterative optimization and upgrading of the risk model. According to the fracturing fracture length and the stress disturbance range, the reservoir reconstruction effect and the prevention of fracturing and casing deformation accidents are taken into account to realize the dynamic optimization and adjustment of the fracturing parameters, mainly including (perforation cluster spacing optimization, fracturing scale / displacement optimization, and well spacing). Figure 7 is the fracturing fracture length under the conditions of single-stage fracturing 4 clusters and 8 clusters. Figure 8 is the fracturing fracture propagation result of the natural fracture development type and the fracture development type section. Figure 9 is the A platform fracturing parameter optimization and control chart. This chart is developed on the basis of considering the A platform fine fracture model and stress model to carry out platform fracturing parameter optimization design while preventing and controlling casing deformation.

[0048] The method of the present patent technology establishes a full three-dimensional real fracture model (with dip angle and strike) on the basis of microseismic fracture characterization technology, and based on this, a calculation method integrating a three-dimensional activation risk model, a three-dimensional fracture slip amount model and a three-dimensional casing stress deformation is constructed under a real three-dimensional fracture, so as to realize rapid solution of real casing deformation. Furthermore, the method combines fracturing construction and casing deformation risk prediction organically, completes iterative optimization and upgrading of the risk model, realizes dynamic optimization and adjustment of fracturing parameters, and realizes subsequent fracturing construction decision-making.

[0049] The technical features of the above embodiments can be combined in any manner. In order to make the description simple, all possible combinations of the technical features in the above embodiments are not described, however, as long as the combinations of the technical features do not exist contradictory, they should be considered as the range disclosed by the present application.

Claims

1. A method for calculating deformation of shale reservoir fracturing casing, characterized in that: include: Acquire microseismic data, imaging logging data, core observation data, and drilling and logging data generated by shale reservoirs during fracturing operations; Construct a fracture network model based on microseismic data, imaging logging data, and core observation data; Constructing a geomechanical model based on the fracture network model, well logging data, and drilling and logging data; the geomechanical model includes the magnitude and direction of the three principal stresses of each fracture surface; The direction cosine matrix is ​​calculated based on the angle between the normal line of each fracture surface and the direction of the three principal stresses. The activation risk value of each fracture surface is calculated based on the magnitude of the three principal stresses and the elements in the direction cosine matrix. The fracture surface corresponding to the activation risk value greater than the preset threshold is regarded as a high-risk fracture surface. The high-risk fracture surface is meshed into three-dimensional triangular surface mesh elements, and a system of linear algebraic equations with the shear displacement continuity invariants of the three-dimensional triangular surface mesh elements as unknowns is established to obtain the three-dimensional slip distribution of the high-risk fracture surface. Taking the three-dimensional slip distribution as the displacement boundary condition, the linear hardening elastoplastic constitutive equation is introduced, and the casing mesh model is constructed using eight-node hexahedral elements. The nonlinear deformation of the casing is obtained by iterative calculation of the casing mesh model, and a deformation cloud diagram of the casing inner diameter is drawn.

2. The method for calculating deformation of shale reservoir fracturing casing according to claim 1, characterized in that: The method also includes iteratively updating the fracture network model and the geomechanical model based on the newly added microseismic data, and dynamically optimizing the perforation cluster spacing, fracturing scale / displacement, and well spacing parameters in combination with the fracturing fracture length, stress perturbation range, and casing change risk prediction results.

3. The method for calculating deformation of shale reservoir fracturing casing according to claim 1, characterized in that: The construction of the fracture network model based on the microseismic data, imaging logging data and core observation data specifically includes: Determine the main fracture trend surface and occurrence information of the fracture in the horizontal section of the fracturing well based on microseismic data; the occurrence information includes: fracture surface inclination, fracture surface tendency and fracture surface area; Based on microseismic data, imaging logging data and core observation data, an initial fracture network model of the work area was established; Determining vertical wellbore fracture information based on the imaging logging data and the core observation data; The main fracture trend surface and occurrence 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 of the work area to obtain the fracture network model.

4. The method for calculating deformation of shale reservoir fracturing casing according to claim 1, characterized in that: The constructing of a geomechanical model based on the fracture network model, well logging data and drilling and logging data specifically includes: Based on the well-seismic fusion technology, the structural surface distribution trend of 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 logging data. Based on the data obtained from indoor 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 property model of the work area. The stress conditions of the model boundary were determined by trial and error method, and the fracture network model, the structural model of the work area and the work area attribute model were superimposed to obtain the geomechanical model of the work area.

5. The method for calculating deformation of shale reservoir fracturing casing according to 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 three-dimensional stress specifically includes: According to the angle between the normal line of the fracture surface and the direction of the maximum horizontal principal stress and the angle between the normal line of the fracture surface and the direction of the vertical principal stress, the angle between the normal line of the fracture surface and the direction of the minimum horizontal principal stress is calculated. The specific formula is: ； Where, is the angle between the normal line of the fracture surface and the direction of the minimum horizontal principal stress, is the angle between the normal line of the fracture surface and the direction of the maximum horizontal principal stress, is the angle between the normal line of the fracture surface and the direction of the vertical principal stress; According to the angle between the normal line of the fracture surface and the direction of the minimum horizontal principal stress and the angle between the normal line of the fracture surface and the direction of the vertical principal stress, the direction cosine matrix is ​​constructed. The formula is: 。 6. The method for calculating deformation of shale reservoir fracturing casing according to claim 5, characterized in that: The calculation formula for the activation risk of the fracture surface is: ； in, ； Where, is the activation risk of the fracture surface, is the direction cosine matrix Elements, , p is the pore pressure inside the fracture, 、 and are the maximum horizontal principal stress, vertical principal stress, and minimum horizontal principal stress, respectively.

7. The method for calculating deformation of shale reservoir fracturing casing according to claim 1, characterized in that: The method of establishing a linear algebraic equation system with the shear displacement continuous invariant of the three-dimensional triangular surface mesh unit as the unknown number and solving it to obtain the three-dimensional slip distribution of the high-risk fracture surface specifically includes: placing a plurality of shear displacement discontinuity elements on the high-risk fracture surface; Based on the displacement discontinuity elements and the shear and normal boundary stresses, a linear equation system of the shear displacement continuity invariant of each three-dimensional triangular surface mesh element is constructed, which is: ； Where, 、 and It is Shear and normal boundary stresses on each element, subscripts 1, 2, 3 represent the three directions of local coordinates; 、 and It is Discontinuous components of shear and normal displacements for each element; ( , =1,2,3) are the boundary influence coefficients.

8. The method for calculating deformation of shale reservoir fracturing casing according to claim 1, characterized in that: The three-dimensional slip distribution is used as the displacement boundary condition, the linear hardening elastoplastic constitutive equation is introduced, and the casing mesh model is constructed using eight-node hexahedral elements. The nonlinear deformation of the casing is obtained by iterative calculation of the casing mesh model, and a deformation cloud diagram of the casing inner diameter is drawn, which specifically includes: The initial casing size model is divided into casing mesh models by eight-node hexahedral elements; The fracture shear slip is used as the displacement boundary condition of the casing model, and the linear hardening model is introduced to characterize the elastic-plastic characteristics of the casing. The stiffness matrix and Gaussian numerical integral are used to calculate the stress of the casing unit node, and the node stress and strain are calculated. The three-dimensional cloud diagram of the casing deformation is drawn, and the deformation result curve is extracted to obtain the casing inner diameter deformation.

9. The method for calculating deformation of shale reservoir fracturing casing according to claim 2, characterized in that: The dynamic optimization of perforation cluster spacing, fracturing scale / displacement, and well spacing parameters specifically includes: After each iteration, the multi-point casing deformation cloud map is compared with the measured wellbore logging data to calculate the error. If the error is greater than a set threshold, the perforation cluster spacing, fracturing scale, displacement rate, and well spacing are adjusted, and all actions from updating the fracture network model to calculating the casing deformation 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