A hydraulic fracture trans-layer propagation simulation method based on DDM-perturbation iteration coupling

Abstract

The application discloses a hydraulic fracture translayer extension simulation method based on DDM-perturbation iteration coupling, belongs to the field of oil and gas field development engineering, and aims to solve the problem of insufficient simulation precision of asymmetric translayer extension of a horizontal well hydraulic fracture in stratified non-uniform rock mass. The method is based on a plane strain assumption, constructs a stratified rock mass model containing stratified fracture toughness, establishes a mechanical equation through a displacement discontinuity method (DDM), introduces a perturbation iteration method to dynamically adjust flow distribution, automatically adapts to new rock layer fracture toughness to update and expand criteria during translayer, and outputs vertical extension form and parameters. The method can be used for simulating translayer behavior of a horizontal well vertical fracture, improving calculation precision and efficiency, and providing a reliable tool for fracturing design of shale gas reservoirs and the like.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas field development engineering, specifically relating to the calculation method, system and medium for the cross-layer propagation of vertical hydraulic fractures in horizontal wells in layered non-uniform rock masses. Background Technology

[0002] Hydraulic fracturing is a core technology in unconventional energy development, and the precise control of fracture propagation directly determines oil and gas recovery rates. Unconventional reservoirs generally exhibit layered heterogeneity, with abrupt differences in mechanical parameters such as geostress, elastic modulus, and fracture toughness among the overlying strata, intermediate reservoir (the layer where the initial fracture is located), and underlying strata. This leads to an asymmetric characteristic of "progressing on one side and being blocked on the other" when hydraulic fractures propagate across layers, posing a challenge to predicting fracturing effectiveness.

[0003] Crack propagation is the stress intensity factor at the crack tip. With rock fracture toughness The equilibrium process: when As the cracks expand, the layered rock mass... A step change will disrupt this balance, causing differences in the propagation rate of fracture branches; these differences further lead to changes in the pressure gradient within the fracture and an imbalance in the distribution of fracturing fluid flow, forming a positive feedback loop of "propagation difference - flow redistribution - asymmetric exacerbation", which intensifies the complexity of propagation.

[0004] Current mainstream methods for calculating hydraulic fracture propagation have significant technical limitations and are insufficient to meet the simulation requirements of layered, non-homogeneous rock masses: Classical models such as PKN and KGD are suitable for simulating fracture evolution in homogeneous reservoirs; the Extended Finite Element Method (XFEM) characterizes the singularity of fracture tips through enrichment functions; the Discrete Element Method (DEM) has advantages in simulating fracture initiation and penetration; and the Displacement Discontinuity Method (DDM) is highly efficient in calculating multi-fracture interactions.

[0005] However, in scenarios of asymmetric cross-layer propagation of hydraulic fractures, the above method has the following technical drawbacks: 1. The propagation criterion cannot be dynamically adapted: interlayer fracture toughness 1. A step change occurs, leading to a sudden change in the threshold of the expansion criterion; 2. Asymmetric fluid-structure coupling is difficult to capture: the difference in crack branch propagation velocity triggers a positive feedback loop of "velocity difference - pressure gradient - flow distribution", which existing methods struggle to dynamically capture; 3. Insufficient characterization of energy conversion laws: the energy conversion rate at the moment of cross-layer penetration is due to... The differences are significant, and existing models struggle to quantitatively characterize the laws governing energy transfer and dissipation. These limitations make it difficult for current methods to accurately simulate the asymmetric cross-layer propagation of hydraulic fractures in layered rock masses.

[0006] To address the above problems, this invention provides a simulation method based on the plane strain assumption: constructing a layered... The layered rock mass model is used to establish mechanical equations through DDM, and a perturbation iteration method is introduced to dynamically adjust the flow distribution, automatically adapting to new rock layers when crossing layers. The updated extended criteria ultimately achieve fully coupled simulation of fluid-solid-fracture, improving computational accuracy and efficiency. Summary of the Invention

[0007] In view of this, the purpose of this invention is to provide a simulation method for asymmetric propagation of hydraulic fracturing in layered rock masses, which can simulate the asymmetric cross-layer propagation of fractures in layered rock masses, and provide a reliable numerical simulation tool for the design and construction optimization of hydraulic fracturing projects.

[0008] This invention provides a simulation method for asymmetric propagation of hydraulic fracturing in layered rock masses, the simulation method comprising:

[0009] Determine the parameters of the layered reservoir to be simulated, including the mechanical parameters of the layered rock mass, the geometric parameters of the rock strata, the in-situ stress parameters, the fracturing fluid parameters, the construction parameters, and the auxiliary parameters for model calculation;

[0010] A two-dimensional Cartesian coordinate system is constructed based on the plane strain assumption. The crack is discretized into uniform grid elements along the height direction and classified into crack tip elements, crack boundary elements at the crack wellbore, and crack middle elements.

[0011] Constructed using the displacement discontinuity method Stiffness matrix Derivation of stiffness coefficient The expression establishes the vector of residual pressure within the seam. Discontinuity displacement vector of the crack The relational expression, the elasticity relation equation is: ;

[0012] Based on the fracturing fluid flow control equation and elastic relationship equation, a structure is constructed using the fracture discontinuity displacement vector. For unknown variables The nonlinear equation system was solved using Newton's iterative method, and a perturbation iterative method was introduced to dynamically adjust the flow distribution coefficients of the left and right branch cracks. ;

[0013] Criteria used in the symmetric extension phase The asymmetric expansion phase employs a criterion. and Determine whether the crack has propagated and update the crack geometry parameters;

[0014] Output fracture geometry parameters, fracture mechanical parameters, and fracturing fluid flow parameters;

[0015] Among them, the layered reservoir is divided into the overlying rock layer, namely rock layer 1, the intermediate rock layer, namely rock layer 2, which is the layer where the initial fracture is located, and the underlying rock layer, namely rock layer 3. DDM-perturbation iterative coupling refers to characterizing the elastic relationship between fracture width and pressure through the displacement discontinuity method, and combining it with the perturbation iterative method to realize the dynamic distribution of fracturing fluid flow during asymmetric propagation, and finally achieve the full coupling simulation of fluid-solid-fracture.

[0016] Optionally, the parameters to be determined for the layered reservoir to be simulated specifically include:

[0017] Mechanical parameters of layered rock masses, including the elastic modulus of each rock layer. Poisson's ratio shear modulus and Type I fracture toughness, the Type I fracture toughness of rock layer No. 1 is: The type I fracture toughness of rock stratum No. 2 is The type I fracture toughness of rock stratum No. 3 is Rock strata geometric parameters: the thickness of rock stratum 1 is infinite, and the thickness of rock stratum 2 is... The thickness of rock layer 3 is infinite, and the initial fracture height is... ,Right now The total height of the crack is Initial height of the left branch Initial height of right branch The initial tips are all located in rock layer No. 2; geostress parameters are used to determine the minimum horizontal principal stress of the reservoir. Used to calculate the residual pressure inside the joint. ,in , The fracturing fluid pressure within the fracture; fracturing fluid parameters: the fracturing fluid is a Newtonian fluid with a viscosity of [missing value]. ,definition Used to simplify the flow control equations; uniform hydrostatic pressure within the gap at the initial moment. , hour Construction parameters: constant total pumping rate in the fracturing well. Left branch crack pumping flow rate Right branch fracture pump flow rate ,and Auxiliary parameters for model calculation, grid cell length Unit half length Time step Numerical iteration allowable error , wellbore radius .

[0018] Optionally, the crack mesh generation specifically includes:

[0019] Construct a two-dimensional Cartesian coordinate system, rotate the vertical hydraulic fracture propagation model 90° clockwise, with the fractured well as the origin, and the fracture height direction as... Axis, left branch along Extending in the negative direction of the axis, the right branch along Extending in the positive direction of the axis, the crack width direction is axis;

[0020] The crack is discretized into uniform mesh elements along its height direction, using an element length of... The left and right branch cracks are discretized using a uniform mesh, with a total number of elements. The left branch unit is numbered as The right branch unit is numbered as follows and Each unit is represented by its center point, and the vertical axis is... The x-axis is or Each unit stores information including unit number, location coordinates, and crack width. and In-fracturing hydraulic pressure ;

[0021] Crack tip unit, number For the right branch tip unit, For the left-hand tip unit, the tip velocity must be satisfied. and The tip position is and Tip seam width At the same time, it satisfies the extended criterion of step S5;

[0022] Fracture boundary unit at fractured wellbore, numbered For the right branch wellbore boundary unit, For the left branch wellbore boundary element, the wellbore flow conservation condition must be satisfied. ,and The location is , The location is , The width of the crack at the wellbore;

[0023] Crack intermediate unit, number For the right branch middle unit, As the middle element of the left branch, it only needs to satisfy the fracturing fluid flow control equation and the displacement discontinuity method elastic relationship equation, without any additional special boundary constraints.

[0024] Optionally, the establishment of the elastic relationship matrix specifically includes:

[0025] Build Stiffness matrix This includes matrix dimensions and core parameters; the stiffness matrix has the following dimensions: Wei, that is Square array Consistent with the total number of elements in step S2, the parameters required for construction include the rock layer shear modulus. Poisson's ratio Unit half length Unit position coordinates That is, the first The horizontal coordinate of each unit That is, the first The x-coordinate of each unit;

[0026] Derivation of stiffness coefficient expression, The physical meaning of is the first When the unit generates a unit discontinuous displacement, it is necessary to wait until the [number]th element generates the [number]th [ The load applied to each element is expressed as: The stiffness coefficient calculation between the left and right support elements only requires replacing the position coordinates of the corresponding elements. or ;

[0027] Establish the residual pressure vector within the seam Discontinuity displacement vector of the crack The relational expression, the elasticity relation equation is: , Let be the vector of the remaining pressure within the crack, the first... The elements are , Let be the displacement vector of the crack discontinuity, the th element ;

[0028] The matrix expansion is as follows: ;

[0029] Stiffness matrix The boundary conditions of the three types of units in step S2 need to be matched. When the crack propagates across layers, due to the different rock layers... Unification, solely through updating fracture toughness The extended criterion needs to be adjusted, therefore no change is required. Structure and The expression.

[0030] Optionally, the coupled solution of the nonlinear equation system specifically includes:

[0031] Based on the fracturing fluid flow control equation and elastic relationship equation, a structure is constructed using the fracture discontinuity displacement vector. For unknown variables A system of nonlinear equations, with As the only unknown variable, , build univariate closed nonlinear equation system The one-dimensional flow control equation for fracturing fluid is neglected. Directional flow The equation is Nonlinear equations for each unit From the above governing equations and , Coupling elimination The following elements were subsequently obtained: crack tip element, right branch tip element. Select For the difference-separated scatter points, the discretization scheme is as follows:

[0032] ;

[0033] Left branch tip unit The discrete format only requires replacing the cell number. ; Crack boundary unit at the wellbore Right branch wellbore unit introduces flow distribution coefficient The right branch flow rate is The proportion of left branch flow is The discrete format is:

[0034] ;

[0035] Left branch shaft unit Discrete format replacement unit number is The flow item is ; Crack middle unit, right branch middle unit Using standard central difference, the discrete format is as follows:

[0036] ;

[0037] Left branch middle unit Discrete format replacement unit number is ;

[0038] Newton's iterative method solves the system of equations, using the previous time step. Solution Based on, estimate Initial Iterative Solution at Time , The iteration number; calculate the function vector. Jacobian matrix , for 1 / 2 square matrix, elements Derivation by unit type; convergence judgment, if... and ,but for The numerical solution of the sand is obtained, and the iteration terminates; iterative updates are performed according to... Calculate the new solution and return to step 2; adjust using the perturbation iteration method. The disturbance is introduced, and the asymmetric expansion stage is set. , For small perturbations, the previous moment Estimate the current initial value During the first asymmetric expansion Coupled iteration, Substituting the wellbore element discretization scheme and combining it with Newton's iteration, the actual flow rate of the right branch is calculated according to Poiseuille's principle; convergence control, if , To allow for error in flow iteration, take , If the result is accurate, then it is the exact value; otherwise, iterate again.

[0039] Optionally, the application of crack propagation criteria specifically includes:

[0040] In the symmetrical extension stage, the cracks were confined to layer 2, with both the left and right tips located within layer 2; the criterion was that... Traffic allocation is as follows: ,Right now ;

[0041] During the asymmetric propagation stage, the fracture penetrates the strata, entering either stratum 1 (left branch) or stratum 3 (right branch), and Based on the judgment criteria: That is, the left branch penetrates the stratum to the No. 1 rock layer, and That is, the right branch penetrates the stratum to the No. 3 rock layer;

[0042] Optionally, the results can be output dynamically in real time, including:

[0043] Real-time dynamic output, every time step Output after convergence: crack geometry parameters, left support height Right support height Total height Crack width in each unit Crack mechanical parameters, residual pressure of each element Residual pressure at the wellbore Stress intensity factor at the tip of the left branch Stress intensity factor at the tip of the right branch Fracturing fluid flow parameters: constant total pumping rate in the fracturing well Left branch crack pumping flow rate Right branch fracture pump flow rate Flow distribution coefficient of left and right branch cracks ;

[0044] This application also provides a computer system for implementing a simulation method for hydraulic fracture propagation across layers based on DDM-perturbation iterative coupling, characterized in that it includes:

[0045] Data input module: Receives layered rock mass mechanical parameters, rock stratum geometric parameters, in-situ stress parameters, fracturing fluid parameters, construction parameters, and model calculation auxiliary parameters required for hydraulic fracture propagation simulation. It supports manual input or import of parameter files in Excel / CSV format and automatically verifies parameter integrity.

[0046] Mesh building module: Performs crack mesh generation operation, automatically generates coordinate system, discrete uniform mesh element, and three types of elements, and outputs an element information table, including number, location coordinates, and element type;

[0047] Matrix calculation module: Performs the operation of building a flexible relation matrix based on the input. Element coordinates, automatically calculate stiffness matrix All elements Supports exporting in matrix format;

[0048] Nonlinear solution module: Integrates Newton's iteration method and perturbation iteration method, with a built-in Jacobian matrix calculation subroutine, automatically performing equation system solutions and... Adjustments were made to display the iterative convergence process in real time.

[0049] Extended Criterion Module: Performs crack extension criterion judgment operation, automatically identifies the extension stage (i.e., symmetrical or asymmetrical), updates crack height and unit position, and triggers the recording of hindrance-breakthrough energy conversion data.

[0050] Results output module: Performs simulation results output operations to generate time series curves and numerical tables;

[0051] Optionally, the nonlinear solution module also includes a parallel computing unit and supports breakpoint resume, meaning that the computation can be restarted from the most recently converged time step after an unexpected interruption.

[0052] This application also provides a computer-readable storage medium storing a computer program, characterized in that, when the program is executed by a processor, it sequentially implements a hydraulic fracture penetration propagation simulation method based on DDM-disturbance iterative coupling, including parameter reading and verification, mesh construction, stiffness matrix calculation, nonlinear equation solving (i.e., Newton iteration combined with disturbance iteration), extended criterion judgment, result output and format conversion, and the program has a built-in parameter anomaly handling subroutine, which prompts for supplementation when parameters are missing and provides reference values ​​when parameters exceed a reasonable range.

[0053] Thus, the fully coupled simulation method and apparatus for hydraulic fracturing and cross-layer propagation of layered rock masses provided in this application embodiment includes: acquiring input parameters of the target layered reservoir, wherein the parameters cover the elastic modulus, Poisson's ratio and fracture toughness, geometric thickness and initial fracture height of each rock layer, in-situ stress and fracturing fluid viscosity, total pump flow rate, and auxiliary quantities for grid and time step calculation; constructing a coordinate system based on the plane strain assumption and uniformly discretizing the left and right branch fractures along the height direction, distinguishing the tip unit, wellbore boundary unit and intermediate unit; and establishing the elastic relationship equation between the residual pressure and discontinuous displacement within the fracture based on the displacement discontinuity method. The stiffness matrix elements are calculated; the one-dimensional fracturing fluid flow control equation is coupled with the elastic relationship equation to construct a nonlinear equation system with element displacement as the unknown, which is solved using Newton's iteration method. Simultaneously, a perturbation iteration method is introduced to dynamically adjust the flow distribution coefficients of the left and right branches. Based on the stress intensity factor criterion, the symmetric stage is taken. Test each asymmetric stage Based on this, it determines whether the fracture has penetrated the layer and updates the fracture geometry in real time; it outputs the fracture height and width distribution, residual pressure and wellbore pressure, and flow rates of the left and right branches in the time step sequence. Results: The above methods enable the coupled solution of the fluid-solid-fracture process and the prediction of cross-layer behavior, providing a quantitative basis for fracturing parameter optimization and risk control. Attached Figure Description

[0054] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0055] Figure 1 This is a schematic diagram of the hydraulic fracture propagation simulation algorithm based on DDM-perturbation iterative coupling provided in an embodiment of the present invention, showing the flow distribution coefficient. The iterative solution steps include initial value estimation, Newton iteration, and convergence judgment process.

[0056] Figure 2 The three-dimensional schematic diagram of the vertical hydraulic fracture provided in the embodiment of the present invention shows the geometry of the vertical hydraulic fracture in the reservoir and the basis for simplifying it into a plane strain problem;

[0057] Figure 3 The following are schematic diagrams of planar hydraulic fracture propagation in reservoirs of different properties provided in the embodiments of the present invention: (a) shows the symmetrical propagation morphology of fractures in a homogeneous reservoir, and (b) shows the asymmetrical propagation morphology of fractures in a layered reservoir.

[0058] Figure 4 This is a comparative schematic diagram of two hydraulic fracture propagation modes provided in the embodiments of the present invention. (a) is the symmetrical fracture propagation mode in a homogeneous reservoir, and (b) is the asymmetrical fracture propagation mode in a layered heterogeneous reservoir.

[0059] Figure 5 This is a schematic diagram of the one-dimensional fracture grid unit division strategy provided in an embodiment of the present invention, showing the unit classification (tip unit, wellbore boundary unit, intermediate unit) and numbering rules of the fractured well and the upper and lower support fractures; Detailed Implementation

[0060] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. The components of the embodiments of the present invention described and labeled in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0061] Hydraulic fracturing is a core technology in unconventional energy development, and the precise control of fracture propagation directly determines oil and gas recovery rates. Unconventional reservoirs generally exhibit layered heterogeneity, with abrupt differences in mechanical parameters such as geostress, elastic modulus, and fracture toughness among the overlying strata, intermediate reservoir (the layer where the initial fracture is located), and underlying strata. This leads to an asymmetric characteristic of "progressing on one side and being blocked on the other" when hydraulic fractures propagate across layers, posing a challenge to predicting fracturing effectiveness.

[0062] Crack propagation is the stress intensity factor at the crack tip. With rock fracture toughness The equilibrium process: when As the cracks expand, the layered rock mass... A step change will disrupt this balance, causing differences in the propagation rate of fracture branches; these differences further lead to changes in the pressure gradient within the fracture and an imbalance in the distribution of fracturing fluid flow, forming a positive feedback loop of "propagation difference - flow redistribution - asymmetric exacerbation", which intensifies the complexity of propagation.

[0063] Current mainstream methods for calculating hydraulic fracture propagation have significant technical limitations and are insufficient to meet the simulation requirements of layered, non-homogeneous rock masses: Classical models such as PKN and KGD are suitable for simulating fracture evolution in homogeneous reservoirs; the Extended Finite Element Method (XFEM) characterizes the singularity of fracture tips through enrichment functions; the Discrete Element Method (DEM) has advantages in simulating fracture initiation and penetration; and the Displacement Discontinuity Method (DDM) is highly efficient in calculating multi-fracture interactions.

[0064] However, in scenarios of asymmetric cross-layer propagation of hydraulic fractures, the above method has the following technical drawbacks: 1. The propagation criterion cannot be dynamically adapted: interlayer fracture toughness 1. A step change occurs, leading to a sudden change in the threshold of the expansion criterion; 2. Asymmetric fluid-structure coupling is difficult to capture: the difference in crack branch propagation velocity triggers a positive feedback loop of "velocity difference - pressure gradient - flow distribution", which existing methods struggle to dynamically capture; 3. Insufficient characterization of energy conversion laws: the energy conversion rate at the moment of cross-layer penetration is due to... The differences are significant, and existing models struggle to quantitatively characterize the laws governing energy transfer and dissipation. These limitations make it difficult for current methods to accurately simulate the asymmetric cross-layer propagation of hydraulic fractures in layered rock masses.

[0065] Based on this, the embodiments of this application provide a simulation method for asymmetric propagation of hydraulic fracturing in layered rock masses, which can simulate the asymmetric cross-layer propagation of fractures in layered rock masses. At the same time, the accompanying computer system and storage medium ensure the implementation of the method, thereby providing a reliable numerical simulation tool for the design and construction optimization of hydraulic fracturing engineering.

[0066] Please see Figure 1 , Figure 1 A schematic diagram of a hydraulic fracture penetration propagation simulation algorithm based on DDM-perturbation iterative coupling provided in an exemplary embodiment of this application is shown.

[0067] like Figure 1 As shown in the exemplary embodiment of this application, the hydraulic fracture propagation simulation method based on DDM-perturbation iterative coupling includes the following steps:

[0068] S101. Determine the parameters of the layered reservoir to be simulated, including the mechanical parameters of the layered rock mass, the geometric parameters of the rock strata, the in-situ stress parameters, the fracturing fluid parameters, the construction parameters, and the auxiliary parameters for model calculation.

[0069] Specifically, the layered rock mass parameters include the elastic modulus of each rock layer. Poisson's ratio shear modulus and Type I fracture toughness The type I fracture toughness of rock stratum No. 1 is The type I fracture toughness of rock stratum No. 2 is The type I fracture toughness of rock stratum No. 3 is ;

[0070] The geometric parameters of the rock strata are as follows: the thickness of rock stratum 1 is infinite, and the thickness of rock stratum 2 is... The thickness of rock layer 3 is infinite, and the initial fracture height is... ,Right now The total height of the crack is Initial height of the left branch Initial height of right branch The initial tips are all located in rock layer No. 2;

[0071] The aforementioned geostress parameters, specifically the minimum horizontal principal stress of the reservoir, are used to calculate the residual pressure within the fracture. , For internal fracturing, hydraulic pressure is strong;

[0072] The fracturing fluid parameters are specified; the fracturing fluid is a Newtonian fluid with a viscosity of [missing value]. ,definition Used to simplify the flow control equations; uniform hydrostatic pressure within the gap at the initial moment. , hour ;

[0073] The construction parameters include a constant total pumping rate in the fracturing well. Left branch crack pumping flow rate Right branch crack pump flow rate ,and ;

[0074] The model calculates auxiliary parameters, including the grid cell length. Unit half length Time step Numerical iteration allowable error , wellbore radius .

[0075] Optionally, the parameters can be adjusted according to actual reservoir conditions, such as elastic modulus. The value is 30 GPa, and the Poisson's ratio is... The initial crack height is set to 0.2. The value is 1m, and the pumping rate is... Value 1×10 -4 m 2 / s, fracturing fluid viscosity Takes a value of 1×10 - 3 Pa·s, initial residual pressure The value is 0.5 MPa, and the element length is... The value is 0.01m, and the time step is... Using a value of 0.01s, calculate the allowable error. Value 1×10 -9 .

[0076] S102. Crack mesh generation: Based on the plane strain assumption, a coordinate system is constructed, and the crack is discretized into uniform mesh elements along the height direction, and classified into crack tip elements, crack boundary elements at the wellbore, and crack middle elements.

[0077] Specifically, the coordinate system construction includes rotating the vertical hydraulic fracture propagation model counterclockwise by 90°, with the fractured well as the origin and the fracture height direction as... Axis, left branch along Extending in the negative direction of the axis, the right branch along Extending in the positive direction of the axis, the crack width direction is The axes form a plane strain rectangular coordinate system;

[0078] The grid cells are discrete, including those with a cell length of... The left and right branch cracks are discretized using a uniform mesh, with a total number of elements. The left branch unit is numbered as The right branch unit is numbered as follows and Each unit is represented by its center point, with the vertical coordinate being 0 and the horizontal coordinate being... or Each unit stores information including unit number, location coordinates, and crack width. and In-fracturing hydraulic pressure ;

[0079] The crack tip unit, numbered For the right branch tip unit, For the left-hand tip unit, the tip velocity must be satisfied. and The tip position is and , tip seam width ;

[0080] The fracture boundary unit at the wellbore is numbered... For the right branch wellbore boundary unit, For the left branch wellbore boundary element, the wellbore flow conservation condition must be satisfied. ,and The location is , The location is , The width of the crack at the wellbore;

[0081] The intermediate unit of the crack is numbered For the right branch middle unit, As the middle element of the left branch, it only needs to satisfy the fracturing fluid flow control equation and the displacement discontinuity method elastic relationship equation, without any additional special boundary constraints.

[0082] Optionally, the meshing ensures that the crack tip elements satisfy strong constraint conditions. and Flow distribution coefficient embedded in wellbore boundary unit The intermediate units employ standard central difference to improve computational accuracy and stability.

[0083] S103. Establishment of the elastic relationship matrix using the displacement discontinuity method. Stiffness matrix Derivation of stiffness coefficient The expression establishes the vector of residual pressure within the seam. Discontinuity displacement vector of the crack Elasticity relation equation ;

[0084] Specifically, matrix dimensions and core parameters, stiffness matrix. for Square array Consistent with the total number of elements in step S2, the parameters required for construction include the rock layer shear modulus. Poisson's ratio Unit half length Unit position coordinates That is, the first The horizontal coordinate of each unit That is, the first The x-coordinate of each unit;

[0085] The matrix dimensions and core parameters, stiffness matrix for Square array Consistent with the total number of elements in step S2, the parameters required for construction include the rock layer shear modulus. Poisson's ratio Unit half length Unit position coordinates That is, the first The horizontal coordinate of each unit That is, the first The x-coordinate of each unit;

[0086] stiffness coefficient Derivation, The physical meaning of is the first When the unit generates a unit discontinuous displacement, it is necessary to wait until the [number]th element generates the [number]th [discontinuous displacement]. The load applied by each element is expressed as: The stiffness coefficient calculation between the left and right support elements only requires replacing the position coordinates of the corresponding elements. or ;

[0087] The elasticity equation is constructed in quantitative form: ,in for The residual pressure vector within the dimensional joint, the first element , for Discontinuity displacement vector of crack, the first element ;

[0088] The matrix expansion is as follows: ;

[0089] Matrix constraints, stiffness matrix The boundary conditions of the three types of units in step S2 need to be matched. When the crack propagates across layers, due to the different rock layers... Unification, solely through updating fracture toughness Adjust the extended criteria; no changes are required. Structure and The expression.

[0090] Optionally, the stiffness matrix The model is constructed based on linear elasticity theory, neglecting the tangential stress between the fracturing fluid and the fracture wall, to ensure computational stability when the fracture tip exhibits strong singularity.

[0091] S104. Nonlinear coupling solution, based on the fracturing fluid flow control equation and elastic relationship equation, constructs a solution... For unknown variables The nonlinear equation system was solved using Newton's iterative method, and a perturbation iterative method was introduced to dynamically adjust the flow distribution coefficients of the left and right branch cracks. .

[0092] Specifically, the nonlinear equation system is constructed based on the one-dimensional flow control equation of fracturing fluid and the elastic relationship equation of the displacement discontinuity method. As the only unknown variable, , build univariate closed nonlinear equation system ;

[0093] The one-dimensional flow control equation for the fracturing fluid is ignored. Directional flow The equation is ;

[0094] The nonlinear equations of each unit From the above governing equations and , Coupling elimination Later obtained;

[0095] The crack tip unit, the right branch tip unit Select For the difference-separated scatter points, the discretization scheme is as follows: Left branch tip unit The discrete format only requires replacing the cell number. ;

[0096] The crack boundary unit at the wellbore. Right branch wellbore unit introduces flow distribution coefficient The right branch flow rate is The proportion of left branch flow is The discrete format is as follows:

[0097] ;

[0098] The left branch shaft unit Discrete format replacement unit number is The flow item is ;

[0099] The crack middle unit, right branch middle unit Using standard central difference, the discrete scheme is as follows: Left branch middle unit Discrete format replacement unit number is ;

[0100] The Newton iteration method is used to solve the problem, with the previous time step... Solution Based on, estimate Initial Iterative Solution at Time , The iteration number; calculate the function vector. Jacobian matrix , for 1 / 2 matrix, elements Derivation by unit type; convergence judgment, if... and ,but for The numerical solution of the sand is obtained, and the iteration terminates; iterative updates are performed according to... Calculate the new solution and return to step 2;

[0101] The perturbation iteration method is adjusted The disturbance is introduced, and the asymmetric expansion stage is set. , For small perturbations, the previous moment Estimate the current initial value During the first asymmetric expansion Coupled iteration, Substituting the wellbore element discretization scheme and combining it with Newton's iteration, the actual flow rate of the right branch is calculated according to Poiseuille's principle; convergence control, if , To allow for error in flow iteration, If the result is accurate, then it is the exact value; otherwise, iterate again.

[0102] Optionally, the Jacobian matrix Calculated based on element type, such as tip element. Including (C.1)-(C.3), ensuring efficient convergence of high-order nonlinear equations; perturbation iteration constraints. To avoid numerical instability.

[0103] S105. Application of crack propagation criteria: Criterion used in the symmetrical propagation stage. All conditions must be met, and the stress intensity factor must be no less than the fracture toughness of the intermediate rock layer. The criterion used in the asymmetric propagation stage is... and It determines whether the crack has expanded and updates the crack geometry parameters.

[0104] Specifically, during the symmetrical extension stage, the cracks were confined to layer 2, with both the left and right tips located within layer 2; the criterion was that... Traffic allocation is as follows: ,Right now ;

[0105] During the asymmetric propagation stage, the fracture penetrates the stratum 1 (left branch) or stratum 3 (right branch), and Based on the judgment criteria: That is, the left branch penetrates the stratum to the No. 1 rock layer, and That is, the right branch penetrates the stratum to the No. 3 rock layer;

[0106] S6: Output results, including fracture geometry parameters, fracture mechanical parameters, and fracturing fluid flow parameters.

[0107] Optionally, the output content may include real-time dynamic output, for each time step. Output after convergence: crack geometry parameters, left support height Right support height Total height Crack width in each unit Crack mechanical parameters, residual pressure of each element Residual pressure at the wellbore Stress intensity factor at the tip of the left branch Stress intensity factor at the tip of the right branch Fracturing fluid flow parameters: constant total pumping rate in the fracturing well Left branch crack pumping flow rate Right branch crack pump flow rate Flow distribution coefficient of left and right branch cracks ;

[0108] Optionally, the output is presented in a visual format, including crack height. , Evolution curve over time, flow distribution Dynamic curves, wellbore pressure and stress intensity factor , Trends in change; supports exporting to tables or charts for easy analysis of energy conversion and asymmetric expansion mechanisms.

[0109] Please see Figure 2 , Figure 3 , Figure 4 , Figure 5 The simulation method provided in this embodiment of the invention includes:

[0110] The computer system is used to implement the above-mentioned simulation method for asymmetric propagation of hydraulic fracturing in layered rock masses, including:

[0111] Data input module: Receives all the above parameters, supports manual input or import of parameter files in Excel / CSV format, and automatically verifies the integrity of the parameters;

[0112] Mesh construction module: Performs the above crack mesh generation, automatically generates three types of elements: coordinate system, discrete uniform mesh element, and classified elements, and outputs an element information table, including element number, location coordinates, and element type;

[0113] Matrix calculation module: Performs the above-mentioned elastic relationship matrix establishment based on the input. Element coordinates, automatically calculate stiffness matrix All elements Supports exporting in matrix format;

[0114] Nonlinear solution module: Integrates the Newton iteration method and the perturbation iteration method mentioned above, with a built-in Jacobian matrix calculation subroutine, automatically performing equation system solutions and... Adjustments are made to display the iterative convergence process in real time;

[0115] Extended Criterion Module: Executes the above extended criteria judgment, automatically identifies the extension stage (i.e., symmetrical or asymmetrical), updates the crack height and unit position, and triggers the recording of the hindrance-breakthrough energy conversion data.

[0116] Result output module: Performs the above result output to generate time series curves and numerical tables;

[0117] Optionally, the nonlinear solution module also includes a parallel computing unit and supports breakpoint resumption, meaning that after an unexpected interruption, the calculation can be restarted from the most recently converged time step.

[0118] The computer-readable storage medium stores a computer program thereon. When the program is executed by the processor, it sequentially implements the steps of the above-mentioned asymmetric propagation simulation method for hydraulic fracturing of layered rock masses, including parameter reading and verification, mesh construction, stiffness matrix calculation, nonlinear equation solving (i.e., Newton iteration combined with perturbation iteration), extended criterion judgment, result output and format conversion. The program also has a built-in parameter anomaly handling subroutine, which prompts for supplementation when parameters are missing and provides reference values ​​when parameters exceed a reasonable range.

[0119] Finally, it should be noted that the above descriptions are merely specific embodiments of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A simulation method for hydraulic fracture propagation across layers based on DDM-perturbation iterative coupling, characterized in that, The simulation method includes the following steps: S1 parameters are determined by identifying the parameters of the layered reservoir to be simulated, including the mechanical parameters of the layered rock mass, the geometric parameters of the rock strata, the geostress parameters, the fracturing fluid parameters, the construction parameters, and the auxiliary parameters for model calculation. S2 fracture mesh generation: Based on the plane strain assumption, a two-dimensional Cartesian coordinate system is constructed. The fracture is discretized into uniform mesh elements along the height direction and classified into fracture tip elements, fracture boundary elements at the fracture wellbore, and fracture middle elements. The S3 elastic relationship matrix is ​​established using the displacement discontinuity method. Stiffness matrix Derivation of stiffness coefficient The expression establishes the vector of residual pressure within the seam. Discontinuity displacement vector of the crack The relational expression, the elasticity relation equation is: ; S4 Nonlinear Equations Coupling Solution: Based on the fracturing fluid flow control equations and elastic relationship equations, a solution is constructed using the fracture discontinuity displacement vector. For unknown variables The nonlinear equation system was solved using Newton's iterative method, and a perturbation iterative method was introduced to dynamically adjust the flow distribution coefficients of the left and right branch cracks. ; The S5 fracture propagation criterion is applied to divide the propagation stage into two phases: symmetrical and asymmetrical, based on the location of the fracture tip within the stratum. The symmetrical propagation stage is characterized by both the left and right fracture tips being located within stratum 2, where the initial fracture was located. The criterion is then applied... The asymmetric propagation stage is characterized by the tip of the left branch of the fracture penetrating the overlying No. 1 rock layer and the tip of the right branch penetrating the underlying No. 3 rock layer. The criteria used are as follows: and In the formula , These are the Type I stress intensity factors at the tips of the left and right branches of the crack, respectively. , , The Type I fracture toughness of rock layers 1, 2 and 3 are respectively used. The propagation criteria of the corresponding stage are used to determine whether the cracks propagate and the crack geometry parameters are updated when the criteria are met. S6 outputs results in real time and dynamically, including fracture geometry parameters, fracture mechanical parameters, and fracturing fluid flow parameters. Among them, the layered reservoir is divided into the overlying rock layer, namely rock layer No. 1, the intermediate rock layer, namely rock layer No. 2, which is the layer where the initial fracture is located, and the underlying rock layer, namely rock layer No.

3. DDM-perturbation iterative coupling refers to characterizing the elastic relationship between fracture width and pressure through the displacement discontinuity method, and combining it with the perturbation iterative method to realize the dynamic distribution of fracturing fluid flow rate during asymmetric propagation, and finally achieve the full coupling simulation of fluid-solid-fracture.

2. The method according to claim 1, characterized in that, Step S1, which involves determining the parameters of the layered reservoir to be simulated, specifically includes: Mechanical parameters of layered rock masses, including the elastic modulus of each rock layer. Poisson's ratio shear modulus and Type I fracture toughness, the Type I fracture toughness of rock layer No. 1 is: The type I fracture toughness of rock stratum No. 2 is The type I fracture toughness of rock stratum No. 3 is ; Rock strata geometric parameters: Rock stratum 1 has an infinite thickness, Rock stratum 2 has a thickness of... The thickness of rock layer 3 is infinite, and the initial fracture height is... ,Right now The total height of the crack is Initial height of the left branch Initial height of right branch The initial tips are all located in rock layer No. 2; In-situ stress parameters are used to determine the minimum horizontal principal stress of the reservoir. Used to calculate the residual pressure inside the joint. ,in , For internal fracturing, hydraulic pressure is strong; Fracturing fluid parameters: The fracturing fluid is a Newtonian fluid with a viscosity of [missing value]. ,definition Used to simplify the flow control equations; uniform hydrostatic pressure within the gap at the initial moment. , hour ; Construction parameters, constant total pumping rate in fracturing wells Left branch crack pumping flow rate Right branch crack pump flow rate ,and ; Model calculation auxiliary parameters, mesh cell length Unit half length Time step Numerical iteration allowable error , wellbore radius .

3. The method according to claim 1, characterized in that, The crack mesh division in step S2 specifically includes: Construct a two-dimensional Cartesian coordinate system, rotate the vertical hydraulic fracture propagation model 90° clockwise, with the fractured well as the origin, and the fracture height direction as... Axis, left branch along Extending in the negative direction of the axis, the right branch along Extending in the positive direction of the axis, the crack width direction is axis; The crack is discretized into uniform mesh elements along its height direction, using an element length of... The left and right branch cracks are discretized using a uniform mesh, with a total number of elements. The left branch unit is numbered as The right branch unit is numbered as follows and Each unit is represented by its center point, and the vertical axis is... The x-axis is or Each unit stores information including unit number, location coordinates, and crack width. and In-fracturing hydraulic pressure ; Crack tip unit, number For the right branch tip unit, For the left-hand tip unit, the tip velocity must be satisfied. and The tip position is and , tip seam width At the same time, it satisfies the extended criterion of step S5; Fracture boundary unit at fractured wellbore, numbered For the right branch wellbore boundary unit, For the left branch wellbore boundary element, the wellbore flow conservation condition must be satisfied. ,and The location is , The location is , The width of the crack at the wellbore; Crack intermediate unit, number For the right branch middle unit, As the middle element of the left branch, it only needs to satisfy the fracturing fluid flow control equation and the displacement discontinuity method elastic relationship equation, without any additional special boundary constraints.

4. The method according to claim 1, characterized in that, Step S3, establishing the elastic relationship matrix, specifically includes: Build Stiffness matrix This includes matrix dimensions and core parameters; the stiffness matrix has the following dimensions: Wei, that is Square array Consistent with the total number of elements in step S2, the parameters required for construction include the rock layer shear modulus. Poisson's ratio Unit half length Unit position coordinates That is, the first The horizontal coordinate of each unit That is, the first The x-coordinate of each unit; Derivation of stiffness coefficient expression, The physical meaning of is the first When the unit generates a unit discontinuous displacement, it is necessary to wait until the [number]th element generates the [number]th [discontinuous displacement]. The load applied by each element is expressed as: The stiffness coefficient calculation between the left and right support elements only requires replacing the position coordinates of the corresponding elements. or ; Establish the residual pressure vector within the seam Discontinuity displacement vector of the crack The relational expression, the elasticity relation equation is: , Let be the vector of the remaining pressure within the crack, the first... The elements are , Let be the displacement vector of the crack discontinuity, the th element ; The matrix expansion is as follows: ; Stiffness matrix The boundary conditions of the three types of units in step S2 need to be matched. When the crack propagates across layers, due to the different rock layers... Unification, solely through updating fracture toughness The extended criterion needs to be adjusted, therefore no change is required. Structure and The expression.

5. The method according to claim 1, characterized in that, Step S4, the coupled solution of the nonlinear equation system, specifically includes: Based on the fracturing fluid flow control equation and elastic relationship equation, a structure is constructed using the fracture discontinuity displacement vector. For unknown variables A system of nonlinear equations, with As the only unknown variable, , build univariate closed nonlinear equation system The one-dimensional flow control equation for fracturing fluid is neglected. Directional flow The equation is Nonlinear equations for each unit From the above governing equations and , Coupling elimination The following elements were subsequently obtained: crack tip element, right branch tip element. Select For the difference-separated scatter points, the discretization scheme is as follows: ； Left branch tip unit The discrete format only requires replacing the cell number. ; Crack boundary unit at the wellbore Right branch wellbore unit introduces flow distribution coefficient The right branch flow rate is The proportion of left branch flow is The discrete format is: ； Left branch shaft unit Discrete format replacement unit number is The flow item is ; Crack middle unit, right branch middle unit Using standard central difference, the discrete format is as follows: ； Left branch middle unit Discrete format replacement unit number is ; Newton's iterative method solves the system of equations, using the previous time step. Solution Based on, estimate Initial Iterative Solution at Time , The iteration number; calculate the function vector. Jacobian matrix , for 1 / 2 matrix, elements Derivation by unit type; convergence judgment, if... and ,but for The numerical solution of the sand is obtained, and the iteration terminates; iterative updates are performed according to... Calculate the new solution and return to step 2; adjust using the perturbation iteration method. The disturbance is introduced, and the asymmetric expansion stage is set. , For small perturbations, the previous moment Estimate the current initial value During the first asymmetric expansion Coupled iteration, Substituting the wellbore element discretization scheme and combining it with Newton's iteration, the actual flow rate of the right branch is calculated according to Poiseuille's principle; convergence control, if , To allow for error in flow iteration, take... , If the result is accurate, then it is the exact value; otherwise, iterate again.

6. The method according to claim 1, characterized in that, The application of the crack propagation criterion in step S5 specifically includes: In the symmetrical extension stage, the cracks extend only within rock layer 2, with both the left and right tips located within rock layer 2; the criterion is that... Traffic allocation is as follows: ,Right now ; During the asymmetric propagation stage, the fracture penetrates the stratum No. 1 (left branch) or stratum No. 3 (right branch), and Based on the judgment criteria: That is, the left branch penetrates the stratum to the No. 1 rock layer, and That is, the right branch penetrates the No. 3 rock layer.

7. The method according to claim 1, characterized in that, The real-time dynamic output of the results in step S6 specifically includes: Real-time dynamic output, every time step Output after convergence: crack geometry parameters, left support height Right support height Total height Crack width in each unit Crack mechanical parameters, residual pressure of each element Residual pressure at the wellbore Stress intensity factor at the tip of the left branch Stress intensity factor at the tip of the right branch Fracturing fluid flow parameters: constant total pumping rate in the fracturing well Left branch crack pumping flow rate Right branch crack pump flow rate Flow distribution coefficient of left and right branch cracks .

8. A computer system for simulating hydraulic fracture propagation across layers based on DDM-perturbation iterative coupling, characterized in that, The system is used to execute the hydraulic fracture propagation simulation method based on DDM-perturbation iterative coupling as described in claim 1. The system includes: Data input module: Receives the mechanical parameters of the layered rock mass, the geometric parameters of the rock strata, the in-situ stress parameters, the fracturing fluid parameters, the construction parameters, and the model calculation auxiliary parameters as described in claim 2. It supports manual input or import of parameter files in Excel / CSV format and automatically verifies the integrity of the parameters. Mesh construction module: Performs crack mesh generation as described in claim 3, automatically generates coordinate system, discrete uniform mesh unit, and three types of units, and outputs a unit information table, including number, location coordinates, and unit type; Matrix calculation module: performs the establishment of the elastic relationship matrix as described in claim 4, based on the input... Element coordinates, automatically calculate stiffness matrix All elements Supports exporting in matrix format; Nonlinear solution module: Integrates the Newton iteration method and perturbation iteration method as described in claim 5, with a built-in Jacobian matrix calculation subroutine, automatically executing the solution of the system of equations and... Adjustments are made to display the iterative convergence process in real time; Extended Criterion Module: Executes the extended criterion judgment as described in claim 6, automatically identifies the extension stage (i.e., symmetrical or asymmetrical), updates the crack height and unit position, and triggers the recording of the hindrance-breakthrough energy conversion data; Result output module: Performs the result output as described in claim 7, generating time series curves and numerical tables.

9. The computer system according to claim 8, characterized in that, The nonlinear solution module also includes a parallel computing unit and supports breakpoint resume, meaning that after an unexpected interruption, the calculation can be restarted from the most recently converged time step.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the program is executed by the processor, it sequentially implements the steps of the hydraulic fracture propagation simulation method based on DDM-disturbance iterative coupling as described in any one of claims 1-7, including parameter reading and verification, mesh construction, stiffness matrix calculation, nonlinear equation system solution (i.e., Newton iteration combined with disturbance iteration), extended criterion judgment, result output and format conversion. The program also has a built-in parameter exception handling subroutine that prompts for supplementation when parameters are missing and provides reference values ​​when parameters exceed a reasonable range.

Images