Geological model-based coal bed gas fracturing three-dimensional well network seam network design method

Abstract

The invention relates to the technical field of geologic model analysis, in particular to a coal bed gas fracturing three-dimensional well network seam network design method based on a geologic model, and the method comprises the steps: firstly constructing a three-dimensional discrete geologic model containing grid node mechanical parameters, and generating a three-dimensional dynamic induced stress potential energy tensor field; performing low-potential-energy gradient optimization calculation by using the field to obtain a shaft space topological coordinate avoiding a high-potential-energy area so as to reduce the well wall instability risk; meanwhile, fracture extension simulation is cut off through an equipotential surface threshold value in the field, and hydraulic fracture geometrical morphology parameters matched with formation energy constraint are accurately determined; and finally, integrating the parameters to generate an optimized three-dimensional well pattern seam network design scheme, and realizing efficient and safe fracturing design based on the rock mass accumulated deformation energy distribution characteristics.

Description

Technical Field

[0001] This invention relates to the field of geological model analysis technology, specifically to a three-dimensional well network and fracture network design method for coalbed methane fracturing based on geological models. Background Technology

[0002] In the exploration and development of unconventional oil and gas resources, especially coalbed methane, hydraulic fracturing technology is often a crucial means to improve reservoir permeability and enhance single-well productivity. Coalbed methane, as a uniquely formed sedimentary rock, often exhibits strong heterogeneity and anisotropy, with a complex internal pore structure. Due to the influence of tectonic movements, the stress distribution and mechanical properties within the formation typically differ significantly in three-dimensional space. Existing well and fracture network design technologies for coalbed methane development largely rely on geological structural interpretation data and conventional rock mechanics parameters. The design focus usually tends to pursue the geometric penetration of the wellbore at geological targets and the volumetric modification range of hydraulic fractures. However, existing design methods, when planning wellbore trajectories, often treat the formation as a relatively homogeneous medium or only consider static stress, and may not adequately quantify the deformation potential energy accumulated by the formation during complex geological histories and the local stress concentration states.

[0003] This situation could lead to the designed wellbore trajectory unknowingly traversing areas of high accumulated deformation energy in the rock mass, thereby increasing the risk of complex engineering accidents such as wellbore instability, collapse, or stuck pipe during drilling. Furthermore, in simulating the extension of hydraulic fractures, existing technologies mostly employ geometric models based on fracture mechanics or simplified extension criteria, lacking sufficient consideration of the physical mechanisms by which the inherent potential energy state of the formation constrains the dynamic extension boundary of fractures. This results in the predicted fracture geometry sometimes failing to accurately reflect the actual limiting effect of formation energy on fracture extension. Therefore, exploring a method that can closely integrate formation stress potential energy distribution with well and fracture network design has significant engineering application value for optimizing coalbed methane development in complex geological environments. Summary of the Invention

[0004] The purpose of this invention is to provide a three-dimensional well network and fracture network design method for coalbed methane fracturing based on geological models, in order to solve the problems mentioned in the background art. Specific technical problems include how to utilize a three-dimensional dynamic induced stress potential energy tensor field constructed based on rock mechanical properties for low-potential-energy gradient path planning and equipotential surface fracture truncation calculation, to solve the technical problems in existing coalbed methane fracturing designs where the wellbore trajectory crosses high-potential-energy instability risk areas due to insufficient consideration of the cumulative deformation energy distribution characteristics of the rock mass, and the insufficient matching degree between the predicted hydraulic fracture morphology and the formation energy state.

[0005] To achieve the above objectives, the present invention aims to provide a three-dimensional well network and fracture network design method for coalbed methane fracturing based on a geological model, specifically including the following method steps:

[0006] S1. Obtain rock mechanical property data for the target coal seam area, specifically including:

[0007] Geophysical logging and core drilling operations were carried out in the target coal seam area. The P-wave transit time, S-wave transit time and logging density of the target coal seam area were collected through geophysical logging. The drilled cores were used to conduct rock mechanics tests in the laboratory to determine the static Young's modulus, static Poisson's ratio and uniaxial compressive strength parameters of the rock.

[0008] Based on the P-wave transit time, S-wave transit time, and logging density, the dynamic mechanical parameters of the rock were calculated. A regression transformation equation was established using the measured static Young's modulus, static Poisson's ratio, and uniaxial compressive strength parameters and the calculated dynamic mechanical parameters of the rock. The regression transformation equation is a univariate linear regression equation established using the measured static mechanical parameters of the rock and the calculated dynamic mechanical parameters of the rock, which is used to correct the dynamic mechanical parameters to static rock mechanical parameters.

[0009] The calculated dynamic rock mechanical parameters are corrected to static rock mechanical parameters using regression transformation equations, resulting in rock mechanical property data of the target coal seam area that are continuously distributed in the depth direction.

[0010] Constructing a three-dimensional discrete geological model that includes the mechanical parameters of the grid nodes, specifically including:

[0011] Using seismic tectonic interpretation data of the target coal seam area, the spatial geometric locations of the coal seam roof, coal seam floor and faults are determined, and a closed three-dimensional geological structural framework is established.

[0012] The three-dimensional geological structure framework is discretized in three-dimensional space to generate a three-dimensional spatial grid system composed of hexahedral elements, and the spatial coordinates of each grid node in the three-dimensional spatial grid system are determined. The three-dimensional spatial grid system is generated by discretizing the three-dimensional geological structure framework using a hexahedral grid partitioning algorithm, and the mechanical parameters of the grid nodes are obtained by mapping the rock mechanical property data to each grid node in the three-dimensional spatial grid system using a kriging interpolation algorithm.

[0013] The Kriging interpolation algorithm is used to map the rock mechanical property data of the target coal seam area to a three-dimensional spatial grid system. Through interpolation calculation, corresponding grid node mechanical parameters are assigned to each grid node, and a three-dimensional discrete geological model containing the grid node mechanical parameters is constructed.

[0014] Step S1 integrates geophysical logging data and laboratory rock mechanics test data, uses regression transformation equations to obtain continuous rock mechanics properties in the depth direction, and maps them to a three-dimensional spatial grid system to construct a three-dimensional discrete geological model. This provides a basic physical model containing accurate spatial geometric information and grid node mechanical parameters for the subsequent construction of a three-dimensional dynamic induced stress potential energy tensor field, ensuring that the energy field calculation is based on real and continuous rock mechanics property data.

[0015] S2. Based on the mechanical parameters of the grid nodes in the three-dimensional discrete geological model, calculate the stress state and cumulative deformation energy of each grid node, and generate a three-dimensional dynamic induced stress potential energy tensor field, specifically including:

[0016] Vertical overlying strata pressure load, horizontal maximum principal stress load, and horizontal minimum principal stress load are applied as boundary conditions to the three-dimensional discrete geological model.

[0017] The finite element numerical simulation algorithm is used to establish the overall stiffness matrix using the mechanical parameters of the mesh nodes and solve the equilibrium equations to calculate the displacement vector of each mesh node.

[0018] The strain tensor is calculated using geometric equations based on the displacement vectors of each grid node, and the triaxial principal stress tensor of each grid node is calculated using physical constitutive equations to determine the stress state of each grid node.

[0019] The cumulative strain energy of each grid node is calculated by substituting the stress state and strain tensor of each grid node into the elastic strain energy density formula.

[0020] The stress state and cumulative deformation energy of each grid node are assigned as tensor properties to the corresponding grid node to generate a three-dimensional dynamic induced stress potential energy tensor field.

[0021] Step S2 applies boundary loads based on a three-dimensional discrete geological model, and uses finite element numerical simulation and elastic strain energy density formula to accurately calculate the stress state and cumulative deformation energy of each grid node, successfully generating a three-dimensional dynamic induced stress potential energy tensor field. This achieves a quantitative characterization of the complex energy distribution inside the stratum rock mass, and solves the problem that it is difficult to accurately assess the potential instability risk area of ​​the stratum due to insufficient consideration of the cumulative deformation energy distribution characteristics of the rock mass.

[0022] S3. Utilizing the three-dimensional dynamic induced stress potential energy tensor field as the path planning space, perform low potential energy gradient optimization calculations to generate wellbore spatial topological coordinates that avoid high potential energy regions. Specifically, this includes:

[0023] The wellhead starting node and the bottom target node are set in the three-dimensional dynamic induced stress potential energy tensor field;

[0024] The cumulative deformation energy values ​​of each grid node in the three-dimensional dynamic induced stress potential energy tensor field are extracted to construct the potential energy scalar field. The spatial rate of change of the cumulative deformation energy values ​​between adjacent grid nodes in the potential energy scalar field is calculated to determine the potential energy gradient vector.

[0025] Grid nodes whose cumulative deformation energy exceeds the preset formation stability threshold are marked as high potential energy regions and set as path planning obstacle points.

[0026] Starting from the wellhead node and ending at the bottom node, the minimum gradient path search algorithm is used to iteratively optimize the grid node set in the non-high potential energy region to select a connected grid node sequence with the minimum cumulative sum of potential energy gradient vector modulus along the way.

[0027] Extract the spatial coordinates of each grid node in the connected grid node sequence in sequence to generate the wellbore spatial topology coordinates that avoid high potential energy regions.

[0028] Step S3 uses the three-dimensional dynamic induced stress potential energy tensor field as the path planning space. By setting high potential energy regions as obstacle points and performing low potential energy gradient optimization calculations, the topological coordinates of the wellbore space that avoid high potential energy regions are selected. This effectively reduces the possibility of the wellbore trajectory crossing high potential energy instability risk regions during the design phase. The energy potential energy guidance mechanism solves the problem of lack of mechanical stability constraints in wellbore path planning.

[0029] S4. Using the wellbore spatial topological coordinates as the fracture initiation point, fracture propagation simulation is performed in a three-dimensional dynamic induced stress potential energy tensor field. The fracture propagation calculation is then truncated using the equipotential surface threshold within the field to determine the geometric parameters of the hydraulic fracture, specifically including:

[0030] Multiple discrete nodes are determined as fracture initiation points from the topological coordinates of the wellbore space;

[0031] In the three-dimensional dynamic induced stress potential energy tensor field, the dominant crack propagation direction perpendicular to the minimum horizontal principal stress is determined based on the stress tensor properties of each grid node; the dominant crack propagation direction is determined by calculating the eigenvalue and eigenvector of the stress tensor at the crack initiation point, and its direction is perpendicular to the eigenvector direction corresponding to the minimum horizontal principal stress.

[0032] A grid node tracing algorithm is used to perform iterative calculations from the crack initiation point to the outer adjacent grid nodes to simulate the crack propagation process in three-dimensional space;

[0033] An equipotential surface threshold is set to characterize the energy boundary at which crack propagation stops, and a spatial equipotential surface with a cumulative deformation energy value equal to the equipotential surface threshold is identified in the three-dimensional dynamic induced stress potential energy tensor field.

[0034] The position of the crack tip is monitored in real time during the iterative calculation of crack propagation simulation, and the crack propagation calculation is immediately cut off when the crack tip touches the spatial equipotential surface.

[0035] Extract the spatial coordinate set of all grid nodes traversed by the crack propagation path, calculate the length, height, and width of the crack based on the spatial coordinate set, and determine the geometric parameters of the hydraulic crack.

[0036] Step S4 simulates fracture propagation in a three-dimensional dynamic induced stress potential energy tensor field and uses an equipotential surface threshold to truncate the fracture extension calculation. By identifying energy equipotential surfaces, the stopping boundary of the fracture is physically defined, thereby determining the geometric parameters of the hydraulic fracture. This method ensures that the prediction of fracture morphology is strictly controlled by the formation energy state, solving the problem of insufficient matching between the predicted hydraulic fracture morphology and the formation energy state in existing designs.

[0037] S5. Based on the wellbore spatial topological coordinates and hydraulic fracture geometry parameters, generate a three-dimensional well network / fracture network design scheme, specifically including:

[0038] A cubic spline interpolation algorithm is used to smoothly fit the topological coordinates of the wellbore space to generate a continuous three-dimensional wellbore trajectory curve;

[0039] Calculate the inclination angle and azimuth angle values ​​of the three-dimensional wellbore trajectory curve at different depth positions;

[0040] Based on the geometric parameters of the hydraulic fractures, the locations of the fracturing segments and perforation clusters are planned on the three-dimensional wellbore trajectory curve.

[0041] A virtual hydraulic fracture mesh model distributed along a three-dimensional wellbore trajectory curve is constructed based on the fracturing segment location, perforation cluster location, and hydraulic fracture geometry parameters.

[0042] The three-dimensional wellbore trajectory curve is spatially combined with the virtual hydraulic fracture mesh model to form a visualized three-dimensional engineering geological entity model;

[0043] The output includes a 3D well network and fracture network design scheme containing well inclination angle values, azimuth angle values, fracturing segment locations, perforation cluster locations, and virtual hydraulic fracture mesh model data.

[0044] Step S5 integrates the optimized wellbore spatial topological coordinates with the hydraulic fracture geometric parameters that conform to energy boundary constraints, constructs a virtual hydraulic fracture mesh model, and generates a visualized three-dimensional well-fracture mesh design scheme. It transforms the field theory analysis results based on rock mechanical properties into specific engineering implementation parameters (such as well inclination angle, azimuth angle, and fracturing segment location), realizing a high-precision well-fracture mesh design that considers the cumulative deformation energy distribution characteristics of the rock mass.

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

[0046] This invention constructs a three-dimensional dynamic induced stress potential energy tensor field based on rock mechanical properties, achieving accurate quantitative characterization of the complex energy distribution within the formation. It automatically plans the wellbore spatial topology coordinates to avoid high-potential-energy risk areas using low-potential-gradient optimization calculations, effectively reducing the risk of wellbore instability and stuck pipe caused by traversing formation stress concentration zones during drilling. Simultaneously, it determines the dynamic propagation boundary of hydraulic fractures based on equipotential surface threshold truncation calculations within the field, ensuring a high degree of matching between the predicted fracture geometry and the actual energy constraint state of the formation. This significantly improves the geological adaptability, engineering safety, and reservoir stimulation efficiency of the three-dimensional well-fracture network design for coalbed methane fracturing. Attached Figure Description

[0047] Figure 1 This is a schematic diagram of the overall method steps of the present invention;

[0048] Figure 2 This is the core flowchart of step S3 of the present invention;

[0049] Figure 3 This is the core flowchart of step S4 of the present invention. Detailed Implementation

[0050] The technical solutions in 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 embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0051] Next, please refer to Figure 1 The purpose of this embodiment is to provide a three-dimensional well network and fracture network design method for coalbed methane fracturing based on a geological model, which includes the following steps:

[0052] S1. Geophysical logging and core drilling operations are carried out in the target coal seam area. P-wave transit time, S-wave transit time, and logging density are collected through geophysical logging. Full-stress-strain rock mechanics testing is conducted on the drilled cores in the laboratory to determine the static Young's modulus, static Poisson's ratio, and uniaxial compressive strength parameters of the rock. Dynamic mechanical parameters of the rock are calculated based on the P-wave transit time, S-wave transit time, and logging density. A regression transformation equation is established using the measured static Young's modulus, static Poisson's ratio, and uniaxial compressive strength parameters and the calculated dynamic mechanical parameters. This regression transformation equation is used to correct the calculated dynamic mechanical parameters to static rock mechanical parameters, resulting in rock mechanical property data of the target coal seam area continuously distributed along the depth direction. The specific formulas for this process are as follows:

[0053] Assume the target coal seam region, as acquired by geophysical well logging, is at a depth of The P-wave time difference at the location is The transverse wave time difference is Density logging values Then calculate the dynamic Poisson's ratio of the rock. With the dynamic Young's modulus of rock The formulas are respectively as well as ,in and Units are uniformly converted to seconds per meter for matching. The International System of Units (SI); subsequently, the static Young's modulus of rock was determined in the laboratory. Static Poisson's ratio and uniaxial compressive strength Establish a univariate linear regression transformation equation: , as well as ,in These represent the static Young's modulus of the rock. Static Poisson's ratio and uniaxial compressive strength The slope coefficient of the regression relationship between them. These represent the intercept constants in the corresponding regression relationships. The calculated dynamic rock mechanical parameters are corrected to static rock mechanical parameters using the above equations, resulting in a sequence of rock mechanical property data continuously distributed along the depth direction. .

[0054] Using seismic tectonic interpretation data of the target coal seam area, the spatial geometric locations of the coal seam roof, floor, and faults are determined, establishing a closed three-dimensional geological structural framework. A hexahedral mesh generation algorithm is used to discretize this framework in three-dimensional space, generating a three-dimensional spatial mesh system composed of hexahedral elements. The spatial coordinates of each mesh node in the three-dimensional spatial mesh system are determined. A kriging interpolation algorithm is used to map the rock mechanical property data of the target coal seam area to the three-dimensional spatial mesh system. Interpolation calculations assign corresponding mechanical parameters to each mesh node, constructing a three-dimensional discrete geological model containing these mechanical parameters. The specific formulas for this process are as follows:

[0055] Define the three-dimensional spatial mesh system generated by discretization as a point set. ,in Representing each grid node, where , and Representing the three-dimensional spatial grid system in , and The total number of mesh nodes along the three coordinate axes; its determined spatial coordinates are... The Kriging interpolation algorithm is used to... Mapped to this system, let the attribute parameters at the grid node to be estimated be... The known attribute values ​​of the valid sample points on the wellbore trajectory are... The interpolation formula is: ,in This represents the total number of valid sample points on the wellbore trajectory used for interpolation calculations, and the weighting coefficients. Unbiased constraints must be satisfied. By solving the system of variational function equations Obtain, among which The variogram matrix is ​​obtained by calculating the pairwise spatial distances between sample points. Indicates the weighting coefficients The vector formed by, and The vector is calculated using valid sample points and the grid nodes to be estimated; thus providing... Each grid node in the array is assigned a corresponding grid node mechanical parameter, denoted as . A three-dimensional discrete geological model containing the mechanical parameters of the grid nodes was constructed.

[0056] The Kriging interpolation algorithm in step S1 is a spatial interpolation method based on the theory of variograms. It is used to predict the attribute values ​​of unsampled points based on the attribute values ​​of known spatial sample points. Its specific implementation process is as follows:

[0057] First, the spatial variability function among valid sample points on the known wellbore trajectory is calculated to quantify the spatial correlation of attribute values ​​as a function of distance, and a variability function model is constructed.

[0058] Next, for each grid node to be estimated in the three-dimensional discrete geological model, a variogram matrix is ​​calculated based on its spatial positional relationship with all known valid sample points and the established variogram model.

[0059] Then, a set of Kriging equations with unbiasedness and optimality constraints is solved. This set of Kriging equations consists of a variogram matrix and a vector composed of the variogram values ​​between known points and points to be estimated, thereby calculating the optimal weight coefficients for each effective sample point.

[0060] Finally, the attribute values ​​of each valid sample point are weighted and summed according to their corresponding optimal weight coefficients to obtain the interpolation estimation result of the grid node to be estimated, thus completing the process of mapping discrete wellbore rock mechanical property data and assigning it to each grid node in the three-dimensional spatial grid system.

[0061] S2. Based on the mechanical parameters of the grid nodes in the three-dimensional discrete geological model, calculate the stress state and cumulative deformation energy of each grid node, and generate a three-dimensional dynamic induced stress potential energy tensor field, specifically including:

[0062] Vertical overburden pressure load, horizontal maximum principal stress load, and horizontal minimum principal stress load are applied as boundary conditions to the three-dimensional discrete geological model. A finite element numerical simulation algorithm is used to establish the overall stiffness matrix using the mechanical parameters of the grid nodes and solve the equilibrium equations. The displacement vectors of each grid node are calculated. Based on the displacement vectors of each grid node, the strain tensor is calculated using geometric equations. Combined with the physical constitutive equations, the triaxial principal stress tensors of each grid node are calculated to determine the stress state of each grid node. The stress state and strain tensor of each grid node are substituted into the elastic strain energy density formula to calculate the cumulative strain energy of each grid node. The stress state and cumulative strain energy of each grid node are assigned as tensor attributes to the corresponding grid nodes, thereby generating a three-dimensional dynamic induced stress potential energy tensor field characterizing the spatial distribution of stress and energy. The specific formulas for this process are as follows:

[0063] Based on a three-dimensional discrete geological model, utilizing the mechanical parameters of the grid nodes In and Constructing the overall stiffness matrix Applying pressure to the overburden layer in the vertical direction Maximum principal stress in the horizontal direction and minimum principal stress boundary load vector By solving the system of linear equilibrium equations The displacement vectors of each grid node were calculated. According to geometric equations Calculate the strain tensor ,in It is a matrix that converts displacement vectors into strain tensors, i.e., the strain-displacement matrix; and combines this with the physical constitutive equations. Calculate the triaxial principal stress tensor for each mesh node. ,in Based on and Define the elastic matrix; then substitute it into the elastic strain energy density formula to calculate the cumulative strain energy of each grid node. ,in The transpose operator represents a matrix or vector; and Assigned as tensor properties to the corresponding mesh nodes to generate a three-dimensional dynamic induced stress potential energy tensor field. .

[0064] The finite element numerical simulation algorithm in step S2 is a numerical calculation method for solving complex engineering mechanics problems. In this method, its specific implementation process is as follows:

[0065] First, based on the geometric information of each hexahedral element in the constructed three-dimensional discrete geological model and the static Young's modulus and static Poisson's ratio in the mechanical parameters of its mesh nodes, the element stiffness matrix of each element is calculated according to the theory of elasticity. This element stiffness matrix characterizes the relationship between the element node force and the node displacement.

[0066] Secondly, according to the spatial position and node number of each element in the overall model, the element stiffness matrices of all elements are assembled into a huge overall stiffness matrix, which represents the mechanical stiffness characteristics of the entire three-dimensional discrete geological model.

[0067] Then, based on the known overlying stratum pressure and horizontal principal stress load conditions applied on the boundary of the three-dimensional discrete geological model, the corresponding nodal load vectors are formed; then, by solving the linear equilibrium equations with the global stiffness matrix as the coefficient matrix, the nodal displacement vector as the unknown quantity, and the nodal load vector as the right-hand side, the displacement vectors of all grid nodes are calculated.

[0068] Finally, based on these displacement solutions, the strain tensor and triaxial principal stress tensor at each grid node are further solved by utilizing the geometric relationship between strain and displacement and the physical constitutive equations of stress and strain, thereby determining its stress state.

[0069] Please see Figure 2 S3. Using the three-dimensional dynamic induced stress potential energy tensor field as the path planning space, perform low potential energy gradient optimization calculations to generate wellbore spatial topological coordinates that avoid high potential energy regions. Specifically, this includes:

[0070] In a three-dimensional dynamic induced stress potential energy tensor field, a wellhead starting node and a bottom-hole target node are set. The cumulative deformation energy values ​​of each grid node in the three-dimensional dynamic induced stress potential energy tensor field are extracted to construct a potential energy scalar field. The spatial rate of change of the cumulative deformation energy values ​​between adjacent grid nodes in the potential energy scalar field is calculated to determine the potential energy gradient vector. Grid nodes whose cumulative deformation energy values ​​exceed a preset formation stability threshold are marked as high potential energy regions and set as path planning obstacle points. Starting from the wellhead starting node and ending at the bottom-hole target node, the minimum gradient path search algorithm is used to iteratively optimize the set of grid nodes in non-high potential energy regions, and select a connected grid node sequence with the minimum cumulative sum of potential energy gradient vector moduli along the path. The spatial coordinates of each grid node in the connected grid node sequence are extracted sequentially to generate the wellbore spatial topology coordinates that avoid high potential energy regions. The specific formula for this process is as follows:

[0071] From the three-dimensional dynamic induced stress potential energy tensor field Extract the cumulative deformation energy values ​​of each grid node. Constructing a potential energy scalar field Calculate the potential energy gradient vector between adjacent grid nodes in the calculation field. ; Set the preset formation stability threshold as The conditions will be met. The grid nodes are defined as the set of obstacle points in the high potential energy region. Starting from the wellhead node Starting from the bottom of the well, the target node Using the minimum gradient path search algorithm as the endpoint, the set of obstacle points in the high potential energy region is searched. Find a connected sequence of grid nodes ,in This represents the total number of grid nodes in the connected grid node sequence; such that the cumulative potential energy gradient vector modulus along the path is summed. Reaching the minimum value, where For spatial step size, This represents the index of the current node when performing a summation calculation in a connected grid node sequence. This refers to the generated wellbore spatial topological coordinates that avoid high-potential-energy regions.

[0072] The minimum gradient path search algorithm in step S3 is an algorithm for searching for the optimal path in a discrete space composed of grid nodes. Its specific implementation process is as follows:

[0073] Starting with the wellhead grid node as the search starting point and the target grid node at the bottom of the well as the search ending point, all grid nodes marked as high-potential-energy obstacle points are excluded from the passable area. In each iteration, from the currently explored node set, the node with the smallest sum of the actual cumulative cost from the starting point to the target node and the estimated heuristic cost from the target node to the end point is selected as the current expansion node. The actual cumulative cost is calculated as the sum of the products of the potential energy gradient vector modulus and the corresponding spatial step size between the adjacent grid nodes traversed at each step on the path from the starting point to the current node. This ensures that the path tends to traverse regions with gentle potential energy changes. By iteratively performing node expansion, neighbor node evaluation, cost update, and path backtracking operations until the target node is found, a connected grid node sequence from the starting point to the end point that effectively avoids high-potential-energy obstacle regions and minimizes the cumulative sum of the potential energy gradient vector modulus along the way is finally selected.

[0074] Please see Figure 3S4. Using the wellbore spatial topological coordinates as the fracture initiation point, fracture propagation simulation is performed in a three-dimensional dynamic induced stress potential energy tensor field. The fracture propagation calculation is then truncated using the equipotential surface threshold within the field to determine the geometric parameters of the hydraulic fracture, specifically including:

[0075] Multiple discrete nodes are identified as fracture initiation points from the wellbore spatial topological coordinates. In a three-dimensional dynamic induced stress potential energy tensor field, the dominant fracture propagation direction perpendicular to the minimum horizontal principal stress is determined based on the stress tensor properties of each grid node. A grid node tracking algorithm is used to iteratively calculate from the fracture initiation point outwards to adjacent grid nodes to simulate the fracture propagation process in three-dimensional space. An equipotential surface threshold is set to characterize the energy boundary at which fracture propagation stops. Spatial equipotential surfaces whose cumulative deformation energy equals the equipotential surface threshold are identified in the three-dimensional dynamic induced stress potential energy tensor field. During the iterative calculation of fracture propagation simulation, the fracture tip position is monitored in real time. When the fracture tip touches the spatial equipotential surface, the fracture propagation calculation is immediately truncated. The spatial coordinate set of all grid nodes traversed by the fracture propagation path is extracted. Based on the spatial coordinate set, the length, height, and width of the fracture are calculated to determine the geometric parameters of the hydraulic fracture. The specific formulas for this process are as follows:

[0076] from Discrete nodes are selected as the crack initiation points. ,exist Extract the stress tensor at that point Then, calculate its eigenvalues ​​and eigenvectors to determine the vector of the dominant crack propagation direction. Perpendicular to the eigenvector direction corresponding to the minimum horizontal principal stress; set the equipotential surface threshold as... The expansion is simulated using a grid node tracking algorithm, with a certain number of iterations. At that time, the position of the crack tip The cumulative deformation energy must satisfy the cutoff criterion: if If the crack propagation path is blocked, the crack will stop expanding; extract the set of spatial coordinates of the mesh nodes traversed by the crack propagation path. ,in This represents the total number of spatial coordinate points of the mesh nodes traversed by the crack propagation path. This represents the discrete fracture node positions that a hydraulic fracture sequentially passes through during its propagation in three-dimensional space, based on the set of spatial coordinates of these grid nodes. Calculate the geometric parameters of the hydraulic fracture, including length. And the height value determined based on the set boundary extrema. and width value ;

[0077] The mesh node tracking algorithm in step S4 is an iterative algorithm used to simulate crack propagation in a discrete mesh system. Its specific implementation process is as follows:

[0078] First, starting from the mesh node that serves as the crack initiation point, the principal stress direction is calculated based on the stress tensor properties of the node in the three-dimensional dynamic induced stress potential energy tensor field, and the direction perpendicular to the minimum horizontal principal stress is initialized as the vector of the crack's dominant extension direction.

[0079] In each iteration step, taking the current crack tip mesh node as the center, among all its adjacent mesh nodes that are not occupied by cracks and are not marked as obstacles, the possible extension direction preference is calculated according to the stress state, and combined with the preset extension criteria (such as the maximum circumferential stress criterion) to select the next mesh node to be extended to, thereby updating the position of the crack tip and the geometric path of the crack.

[0080] During the expansion process, the cumulative deformation energy value at the new tip location is acquired in real time and compared with the preset equipotential surface threshold. Once the value is equal to the equipotential surface threshold, it indicates that the crack tip has touched the spatial equipotential surface representing the energy boundary, and the iterative calculation is terminated immediately, thus simulating the physical process of the crack stopping expansion when it encounters the energy boundary.

[0081] S5. Based on the wellbore spatial topological coordinates and hydraulic fracture geometry parameters, generate a three-dimensional well network / fracture network design scheme, specifically including:

[0082] A cubic spline interpolation algorithm is used to smooth and fit the topological coordinates of the wellbore space to generate a continuous three-dimensional wellbore trajectory curve. The inclination angle and azimuth angle values ​​of the three-dimensional wellbore trajectory curve at different depth positions are calculated. Based on the geometric parameters of the hydraulic fractures, the fracturing segment positions and perforation cluster positions are planned on the three-dimensional wellbore trajectory curve. A virtual hydraulic fracture mesh model distributed along the three-dimensional wellbore trajectory curve is constructed based on the fracturing segment positions, perforation cluster positions, and hydraulic fracture geometric parameters. The three-dimensional wellbore trajectory curve and the virtual hydraulic fracture mesh model are spatially combined to form a visualized three-dimensional engineering geological entity model. The output is a three-dimensional well network fracture network design scheme containing inclination angle values, azimuth angle values, fracturing segment positions, perforation cluster positions, and virtual hydraulic fracture mesh model data.

[0083] The cubic spline interpolation algorithm is used to determine the spatial topological coordinates of the wellbore. Perform smooth fitting to generate a continuous three-dimensional wellbore trajectory curve. ,in For path parameters; calculate the first derivative vector of any point on the trajectory curve. Then, the inclination angle value is obtained. and azimuth values ,in , and Representing the three-dimensional wellbore trajectory curves respectively , and For parameters First derivative components; based on the geometric parameters of hydraulic fractures ,exist The fracturing segment location parameters that satisfy the arc length condition are calculated based on the preset spacing. and perforation cluster position parameters , construct Centered on, the geometric dimensions are determined by Defined virtual hydraulic fracture mesh model The final output consists of data tuples A three-dimensional well network and fracture network design scheme.

[0084] The cubic spline interpolation algorithm in step S5 is a mathematical method for generating smooth curves through a series of discrete data points. Its specific implementation process is as follows:

[0085] The topological coordinates of the wellbore space, i.e., the three-dimensional spatial coordinates of a series of discrete grid nodes connected in sequence, are used as the known type point inputs for the interpolation algorithm. First, a cubic polynomial function is constructed for each curve segment between adjacent type points. This function is defined by the coefficients of the first and second derivatives of the position coordinates with respect to the curve parameters.

[0086] Then, by setting the function values, first derivative values, and second derivative values ​​of adjacent curve segments to be continuous at all internal point values, and applying specific boundary conditions (such as natural boundary conditions), a system of linear equations about all polynomial coefficients is established and solved. After solving this system of linear equations, a smooth three-dimensional spatial curve with continuous second derivatives as a whole is obtained, namely the three-dimensional wellbore trajectory curve. By taking the first derivative of this three-dimensional wellbore trajectory curve, the tangent vector at any point on the curve can be obtained, and then the inclination angle and azimuth angle values ​​corresponding to any depth position on the trajectory can be accurately calculated.

[0087] As can be seen from the above description, the three-dimensional well network and fracture network design method for coalbed methane fracturing based on geological models provided in this embodiment has the following technical effects:

[0088] By establishing a continuous sequence of rock mechanical property data in the depth direction and constructing a three-dimensional discrete geological model containing mechanical parameters of grid nodes, a three-dimensional dynamic induced stress potential energy tensor field characterizing the formation energy distribution was accurately generated. This tensor field was used to perform low-potential-energy gradient optimization calculations, automatically selecting wellbore spatial topological coordinates that avoid high-potential-energy regions, effectively reducing the risk of wellbore instability. Simultaneously, the equipotential surface threshold within the field was used to truncate fracture extension calculations, ensuring that the determined hydraulic fracture geometry parameters are truly constrained by the formation energy boundary. Finally, based on a smooth three-dimensional wellbore trajectory curve and a virtual hydraulic fracture mesh model generated by a cubic spline interpolation algorithm, a visualized three-dimensional well-fracture design scheme was output, significantly improving the geological adaptability, engineering safety, and implementation accuracy of coalbed methane fracturing design.

[0089] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely preferred examples and are not intended to limit the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A three-dimensional well-fracture network design method for coalbed methane fracturing based on a geological model, characterized in that, The methods and steps include the following: S1. Obtain rock mechanical property data of the target coal seam area and construct a three-dimensional discrete geological model containing mechanical parameters of grid nodes; S2. Based on the mechanical parameters of the grid nodes in the three-dimensional discrete geological model, calculate the stress state and cumulative deformation energy of each grid node, and generate a three-dimensional dynamic induced stress potential energy tensor field. S3. Using the three-dimensional dynamic induced stress potential energy tensor field as the path planning space, perform low potential energy gradient optimization calculation to generate wellbore spatial topology coordinates that avoid high potential energy regions. S4. Using the topological coordinates of the wellbore space as the fracture initiation point, the fracture propagation is simulated in the three-dimensional dynamic induced stress potential energy tensor field, and the fracture propagation calculation is truncated by the equipotential surface threshold in the field to determine the geometric parameters of the hydraulic fracture. S5. Generate a three-dimensional well network and fracture network design scheme based on the topological coordinates of the wellbore space and the geometric morphology parameters of the hydraulic fractures.

2. The three-dimensional well network and fracture network design method for coalbed methane fracturing based on a geological model according to claim 1, characterized in that, The process of acquiring the rock mechanical property data specifically includes: Geophysical logging and core drilling operations were carried out in the target coal seam area. The P-wave transit time, S-wave transit time and logging density of the target coal seam area were collected through geophysical logging. The drilled cores were used to conduct rock mechanics tests in the laboratory to determine the static Young's modulus, static Poisson's ratio and uniaxial compressive strength parameters of the rock. Based on the P-wave transit time, S-wave transit time, and logging density, the dynamic mechanical parameters of the rock are calculated, and a regression transformation equation is established using the measured static Young's modulus, static Poisson's ratio, and uniaxial compressive strength parameters and the calculated dynamic mechanical parameters of the rock. The calculated dynamic rock mechanical parameters are corrected to static rock mechanical parameters using the regression transformation equation, resulting in rock mechanical property data of the target coal seam area that is continuously distributed in the depth direction.

3. The three-dimensional well network and fracture network design method for coalbed methane fracturing based on a geological model according to claim 1, characterized in that, The construction process of the three-dimensional discrete geological model specifically includes: Using seismic tectonic interpretation data of the target coal seam area, the spatial geometric locations of the coal seam roof, coal seam floor and faults are determined, and a closed three-dimensional geological structural framework is established. The three-dimensional geological structure framework is discretized in three-dimensional space to generate a three-dimensional spatial grid system composed of hexahedral elements, and the spatial coordinates of each grid node in the three-dimensional spatial grid system are determined. The rock mechanical property data of the target coal seam area are mapped to a three-dimensional spatial grid system using the Kriging interpolation algorithm. The corresponding mechanical parameters of each grid node are assigned to each grid node through interpolation calculation, and a three-dimensional discrete geological model containing the mechanical parameters of the grid nodes is constructed.

4. The three-dimensional well network and fracture network design method for coalbed methane fracturing based on a geological model according to claim 1, characterized in that, The generation process of the three-dimensional dynamically induced stress potential energy tensor field specifically includes: The overlying strata pressure load in the vertical direction, the maximum horizontal principal stress load in the horizontal direction, and the minimum horizontal principal stress load in the horizontal direction are applied as boundary conditions to the three-dimensional discrete geological model. The finite element numerical simulation algorithm is used to establish the overall stiffness matrix using the mechanical parameters of the mesh nodes and solve the equilibrium equations to calculate the displacement vector of each mesh node. The strain tensor is calculated using geometric equations based on the displacement vectors of each grid node, and the triaxial principal stress tensor of each grid node is calculated using physical constitutive equations to determine the stress state of each grid node. The cumulative strain energy of each grid node is calculated by substituting the stress state and strain tensor of each grid node into the elastic strain energy density formula. The stress state and cumulative deformation energy of each grid node are assigned as tensor properties to the corresponding grid node to generate a three-dimensional dynamic induced stress potential energy tensor field.

5. The three-dimensional well network and fracture network design method for coalbed methane fracturing based on a geological model according to claim 1, characterized in that, The process of generating the wellbore spatial topological coordinates specifically includes: In the three-dimensional dynamic induced stress potential energy tensor field, the wellhead starting node and the bottom target node are set; The cumulative deformation energy values ​​of each grid node in the three-dimensional dynamic induced stress potential energy tensor field are extracted to construct the potential energy scalar field. The spatial rate of change of the cumulative deformation energy values ​​between adjacent grid nodes in the potential energy scalar field is calculated to determine the potential energy gradient vector. Grid nodes whose cumulative deformation energy exceeds the preset formation stability threshold are marked as high potential energy regions and set as path planning obstacle points. Starting from the wellhead node and ending at the bottom node, the minimum gradient path search algorithm is used to iteratively optimize the grid node set in the non-high potential energy region to select a connected grid node sequence with the minimum cumulative sum of potential energy gradient vector modulus along the way. Extract the spatial coordinates of each grid node in the connected grid node sequence in sequence to generate the wellbore spatial topology coordinates that avoid high potential energy regions.

6. The three-dimensional well network and fracture network design method for coalbed methane fracturing based on a geological model according to claim 1, characterized in that, The process of determining the geometric parameters of the hydraulic fracture specifically includes: Multiple discrete nodes are determined as fracture initiation points from the wellbore spatial topological coordinates; In the three-dimensional dynamic induced stress potential energy tensor field, the dominant crack propagation direction perpendicular to the minimum horizontal principal stress is determined based on the stress tensor properties of each grid node. A grid node tracing algorithm is used to perform iterative calculations from the crack initiation point to the outer adjacent grid nodes to simulate the crack propagation process in three-dimensional space. An equipotential surface threshold is set to characterize the energy boundary at which crack propagation stops, and a spatial equipotential surface in the three-dimensional dynamic induced stress potential energy tensor field is identified whose cumulative deformation energy value is equal to the equipotential surface threshold. The position of the crack tip is monitored in real time during the iterative calculation of crack propagation simulation, and the crack propagation calculation is immediately cut off when the crack tip touches the spatial equipotential surface. Extract the spatial coordinate set of all grid nodes traversed by the crack propagation path, calculate the length, height, and width of the crack based on the spatial coordinate set, and determine the geometric parameters of the hydraulic crack.

7. The three-dimensional well network and fracture network design method for coalbed methane fracturing based on a geological model according to claim 1, characterized in that, The generation process of the three-dimensional well network and fracture network design scheme specifically includes: A cubic spline interpolation algorithm is used to smooth and fit the spatial topological coordinates of the wellbore to generate a continuous three-dimensional wellbore trajectory curve; Calculate the inclination angle and azimuth angle values ​​of the three-dimensional wellbore trajectory curve at different depth positions; Based on the aforementioned hydraulic fracture geometry parameters, the locations of the fracturing segments and perforation clusters are planned on the three-dimensional wellbore trajectory curve. A virtual hydraulic fracture mesh model is constructed based on the fracturing segment location, the perforation cluster location, and the hydraulic fracture geometry parameters, distributed along the three-dimensional wellbore trajectory curve. The three-dimensional wellbore trajectory curve is spatially combined with the virtual hydraulic fracture mesh model to form a visualized three-dimensional engineering geological entity model; The output includes a three-dimensional well network and fracture network design scheme containing the well inclination angle value, the azimuth angle value, the fracturing segment position, the perforation cluster position, and the virtual hydraulic fracture mesh model data.

8. The three-dimensional well network and fracture network design method for coalbed methane fracturing based on a geological model according to claim 2, characterized in that, The regression transformation equation is a univariate linear regression equation established using the static mechanical parameters of rocks measured in the laboratory and the dynamic mechanical parameters of rocks calculated, which is used to correct the dynamic mechanical parameters to static rock mechanical parameters.

9. The three-dimensional well network and fracture network design method for coalbed methane fracturing based on a geological model according to claim 3, characterized in that, The three-dimensional spatial grid system is generated by discretizing the three-dimensional geological structure framework using a hexahedral grid partitioning algorithm. The mechanical parameters of the grid nodes are obtained by mapping the rock mechanical property data to each grid node in the three-dimensional spatial grid system using a kriging interpolation algorithm.

10. The three-dimensional well network and fracture network design method for coalbed methane fracturing based on a geological model according to claim 6, characterized in that, The dominant direction of crack propagation is determined by calculating the eigenvalue and eigenvector of the stress tensor at the crack initiation point, and its direction is perpendicular to the eigenvector direction corresponding to the minimum horizontal principal stress.

Images