Fracture propagation form prediction method considering gravel characteristics

Abstract

The invention discloses a fracture propagation form prediction method considering gravel characteristics, and the method comprises the steps: building a mesoscale conglomerate hydraulic fracture propagation particle discrete element model considering the real gravel contact behavior based on a Voronoi polygon particle-based algorithm and a fluid flow network; through correction of mechanical parameters and permeability parameters, gravel contact properties and cementation properties of the real conglomerate are matched; the dynamic evolution process of the conglomerate fracturing crack is analyzed, and crack initiation characteristics, a crack asymmetric extension mechanism, a gravel-surrounding extension mechanism and the corresponding relation between the crack width and the crack type are determined; main control factor analysis can be carried out on parameters such as lithology, stress states, gravel properties and engineering parameters, and a real hydraulic fracture-gravel interaction mechanism is disclosed. The device is simple in structure, easy to maintain, low in cost and convenient to popularize. Support is provided for popularization of the stress wall effect in the future, a quantitative characterization means is provided for quantitative evaluation of active stress interference of the reservoir, and theoretical guidance is provided for large-scale three-dimensional development.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas field development, and specifically relates to a method for predicting fracture propagation morphology that takes into account gravel characteristics. Background Technology

[0002] Volumetric fracturing technology for reservoir stimulation can create complex hydraulic fracture networks, increasing reservoir drainage capacity, and is widely used in the development of tight sandstone and conglomerate oil and gas reservoirs. In recent years, with the development and research of unconventional oil and gas resources, the abundant tight sandstone and conglomerate oil and gas reservoirs have received widespread attention and extensive research. However, due to the presence of gravel, the hydraulic fractures are unevenly distributed, and the fracture network usually cannot form sufficient coverage area and complexity. Currently, the understanding of the complex propagation mechanisms of artificial fractures "around gravel" and "through gravel" is not deep enough. Since laboratory experiments cannot directly observe the impact of gravel on fracture propagation behavior during reservoir stimulation, and because laboratory experiments are rare and not reproducible, numerical simulation experiments have become one of the main methods for studying the interference of gravel on fracture propagation in conglomerate reservoirs.

[0003] Current research on hydraulic fracturing numerical simulations is based on the CDEM algorithm: the Discrete Element Method (CDEM) for continuous media couples the Finite Element Method (FEM) and the Discrete Element Method (DEM) to simulate rock fracture. The CDEM algorithm discretizes the rock mass into block elements, each of which can be considered as a single finite element or multiple finite elements. Contact elements are then set between these block elements, as shown in the attached diagram. Figure 2 As shown, each block can be considered as a single finite element or multiple finite elements. The block element internally follows linear elasticity laws and is calculated using the finite element method; while the contact element can be separated using the discrete element method. A two-dimensional contact consists of a normal spring and a tangential spring (see attached diagram). Figure 3 A three-dimensional contact consists of a normal spring and two tangential springs.

[0004] When a material inherently exhibits failure characteristics, transforming into a discontinuous medium, or fails under external loads, thus changing from a continuous to a discontinuous medium, this is represented in CDEM calculations as the fracture of a corresponding spring. When the stress satisfies the maximum tensile stress criterion and the Mohr-Coulomb strength theory, the element contact nodes separate, the contact force becomes zero, and the element itself only deforms without failure. That is, the contact between an element and its adjacent elements can be separated, and the interaction forces between elements are calculated based on the relationship between the force and displacement at the element nodes. The motion and deformation of each block are calculated using the element's force and stiffness matrix. When performing numerical simulations using a mesh model, parameters are assigned to different model meshes based on the matrix mesh and particle mesh division within the model.

[0005] However, current discrete element method (DEM) models for conglomerate hydraulic fracture propagation based on continuous media introduce gravel by assigning different mesh attributes, failing to consider the true linear contact, interlocking, and interfacial slip characteristics between gravel particles. This results in an inability to reproduce the interaction process between hydraulic fractures and the multi-gravel interface, and an inability to accurately understand the microscale evolution mechanism of hydraulic fractures in conglomerate. Patent CN201911010105.4 relates to a numerical simulation method for a complex anisotropic constitutive relation model. This method establishes a discrete element physical model of sandstone and conglomerate by selecting wall elements, generating particles, applying confining pressure, generating a wellbore, and establishing initial equilibrium. It optimizes the fluid-particle coupling mode and contact constitutive relation model to determine the parameters of the sandstone and conglomerate DEM model. The method compares the fracture morphology changes during the fracture penetration process under different confining pressures, gravel strengths, and injection rates. Its advantage lies in the use of a random particle packing pattern, which better reflects the heterogeneity of conglomerate; however, it also has disadvantages: it does not consider the true linear contact, interlocking, and interfacial slip characteristics between gravel particles, and cannot reproduce the interaction process between hydraulic fractures and the multi-gravel interface. Summary of the Invention The purpose of this invention is to provide a method for predicting fracture propagation morphology that takes into account gravel characteristics, in order to solve the problem of how gravel interferes with fracture propagation during hydraulic fracturing of highly heterogeneous sandstone and conglomerate reservoirs.

[0006] The objective of this invention is achieved through the following technical means: a method for predicting crack propagation morphology considering gravel characteristics, comprising the following steps: S1. Based on the Voronoi polygonal grain-based gravel generation algorithm, construct discrete element models of conglomerate particles with different gravel characteristics. Assign smooth joint contact models to the gravel-matrix interface and the gravel-gravel interface, and assign corresponding joint surface contact parameters. Assign linear parallel bonded contact models to the matrix-matrix particle contact and the gravel-gravel particle contact, and assign corresponding bonded parameters. S2, Flowing Time Step Update S2-1, Contact Force Update: Based on the particle contact model, the motion state and force state of the rigid particle system are explicitly iteratively updated. S2-2, Flow channel opening update: Based on the normal contact force on the associated particles, i.e., the force state, update the iterative equation and update the flow channel opening w. p ; S2-3. Obtain the fluid volume change Δ within a single pore domain during a single flow step. V f The fluid volume change Δ is obtained by measuring the injection amount from the source term and the fluid injection amount from the flow channel. V f ; S2-4. Calculate the pore pressure change Δp in each pore of the pore network; S2-5. Determine the steady flow time step and obtain the critical flow time step Δ. t fc ,

[0007] In the formula, m represents the opening of the i-th flow channel connected to the mother pore. Let m be the length of the i-th flow channel connected to the mother pore. K f The bulk modulus of the fluid, Pa. V p For true pore volume, m 2 , This represents the change in pore volume. is the fluid viscosity, mPa·s, and n is the number of particles; Introducing a safety factor α t By replacing the first term on the right-hand side of the above equation, we obtain the steady-state flow time step.

[0008] In the formula, 0 < α t <1, take α t It is 0.8; S3, compared to the steady flow time step Δ t f With stable solid time step Δ t m The smaller of the two values ​​is taken as the final iteration step; S4. Solid time step update: Based on the final iteration time step, repeat step S2 for alternating iterative updates. This involves determining the parameters from the previous solid time step, including particle stress state, particle displacement, particle position, and particle contact force. Finally, update the particle velocity and position, and pore volume V of the current time step. d and flow channel opening w p We determine whether the updated particle state meets the particle equilibrium condition. If it does, it means that the crack propagation is complete. This is used to simulate the hydraulic crack propagation process in a homogeneous medium and obtain numerical simulation results. In S1, the specific steps for constructing discrete element models of conglomerate particles with different gravel characteristics are as follows: Gravel-matrix initial particle generation: Based on the actual gravel size distribution and volume content, a close-packed mass of gravel and matrix initial particles is generated; Gravel-matrix particle contact topology identification: Based on the polygonal loop algorithm, several particles are connected sequentially through the contact head and tail to form a closed loop, which constitutes the basic topological unit; Gravel-matrix polygonal particle topology identification: Connect the centroids of adjacent closed loops to establish the gravel-matrix polygonal particle topology; Gravel-matrix mesh region refilling: Based on the polygonal particle topology, small-sized disk particles are used to refill the model's computational domain. Disk particles located within the gravel mesh region are considered as gravel constituent particles, and disk particles located within the matrix mesh region are considered as matrix constituent particles. At this point, particle contact types can be classified into matrix-matrix particle contact within the matrix, gravel-gravel particle contact within the gravel, gravel-matrix interface particle contact, and gravel-gravel interface particle contact. Smooth joint contact models are assigned to the gravel-matrix interface particle contact and gravel-gravel interface particle contact at the gravel-matrix mesh boundary intersection, and corresponding joint surface contact parameters are assigned; linear parallel bonded contact models are assigned to the matrix-matrix particle contact and gravel-gravel particle contact within the matrix mesh and gravel mesh, and corresponding bonded parameters are assigned.

[0009] In S2-1, the specific method for explicitly iteratively updating the motion and force states of the rigid particle system is as follows: The gravel-matrix particle contact model has been updated to: , ,

[0010] In the formula, F l and The contact forces, in N, are the linear and parallel bonding portions, respectively. The contact torque of the parallel bonding portion is N·m; It is the unit normal vector perpendicular to the plane of the paper and pointing outwards; and N represents the normal and tangential components of the linear partial contact force. and For the normal and tangential components of the contact force in the parallel bonded area, N; M b The bending moment is N·m. For planar normal contact force, N, The tangential contact force is N; The gravel interface contact model has been updated to...

[0011] In the formula, and These are the normal and tangential contact forces on the joint surface, respectively, in N. Let N be the normal contact force in the plane of any element j. Let N be the tangential contact force of any element j.

[0012] In step S2-2, when the contact adhesion is intact and the object is in a compressed state... ; In the formula, w 0 represents the initial opening of the flow channel, in meters (m). F n0 The flow channel opening is w The normal contact force at 0 / 2, in N; when the contact adhesion is intact and under tension, ; When the contact bond breaks ; In the formula, w b The flow channel opening, in meters, corresponds to the instantaneous breakage of the contact bond. m Δ is a dimensionless correction factor for the flow channel opening. w d The change in particle contact gap within the current time step is represented by m, w0, and m, which are variables related to the permeability properties of the pore network, and their values ​​are determined through parameter correction. It is the normal contact force.

[0013]

[0014] In the formula, q p The fluid volumetric flow rate, m 3 / s; w p Let m be the flow channel opening. p 2 and p 1 represents the fluid pressure in pore domain 1 and pore domain 2, respectively, in Pa; μ The viscosity of the fluid is Pa·s; l p Let be the length of the flow channel, in meters (m).

[0015] In S4,

[0016] In the formula, The first pore domain connected to this pore domain i The volumetric flow rate of the injected fluid in each flow channel is m³. 3 / s; The first one connected to this pore domain i Fluid volumetric flow rate in each flow channel, m 3 / s; n This represents the total number of flow channels connected to this pore region; Δ t f1 Let s be the time step during the flow.

[0017] In S2-4, the pore pressure change Δp of each pore in the pore network is specifically as follows:

[0018] In the formula, Δ p The change in pore pressure is expressed in Pa. K f Here, Δ is the fluid bulk modulus, Pa; V p The change in pore volume is m. 2 ; V p For true pore volume, m 2 .

[0019] The beneficial effects of this invention are as follows: Based on the Voronoi polygonal particle-based algorithm and fluid flow network, a microscale discrete element model of hydraulic fracture propagation in conglomerate considering the actual contact behavior of gravel is established; through correction of mechanical and permeability parameters, the gravel contact properties and cementation properties of real conglomerate are matched; the dynamic evolution process of hydraulic fracturing fractures in conglomerate is analyzed, clarifying fracture initiation characteristics, asymmetric fracture propagation mechanism, perigravel propagation mechanism, and the correspondence between fracture width and fracture type; it allows for the analysis of controlling factors such as lithology, stress state, gravel properties, and engineering parameters, revealing the actual interaction mechanism between hydraulic fractures and gravel. This invention has a simple structure, is easy to maintain, has low cost, and is convenient to promote. It can support the future promotion of stress wall effects, provide a quantitative characterization method for evaluating active stress disturbances in reservoirs, provide theoretical guidance for large-scale three-dimensional development, and save on field fracturing test costs. Attached Figure Description

[0020] Figure 1 Flowchart of a crack propagation morphology prediction method that takes into account gravel characteristics; Figure 2 This is a schematic diagram of the block unit of the CDEM algorithm; Figure 3 This is a schematic diagram illustrating the block unit separation of the CDEM algorithm. Figure 4 This is a schematic diagram of particle contact. Figure 5 The crack propagation morphology under different gravel contents; Figure 6 The crack propagation morphology under different gravel particle sizes; The present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Detailed Implementation

[0021]

Example 1

[0022] Based on the discrete element method, a microscale conglomerate fracture propagation unit is constructed using a Voronoi polygon-based real gravel generation algorithm and a fluid flow pore network. The motion and stress states of the rigid particle system within the propagation unit are explicitly iteratively updated according to the particle contact model to simulate the tensile and bending failure processes.

[0023] A microscale hydraulic fracture propagation model for conglomerate is proposed based on the Voronoi polygon-based gravel generation algorithm and a fluid flow pore network. This model comprises multiple conglomerate discrete element models and fluid flow network elements (based on laminar flow equations driven by pressure differences between incompressible fluids in the pore domains at both ends). Using a mineral particle generation algorithm that matches the rock's mineral content and size distribution, the model accurately reflects the topological structure and statistical distribution characteristics of real mineral particles. Discrete element computational domains for conglomerate particles with different gravel contents and sizes are constructed. A pore network model (i.e., joint contact model and linear parallel bonded contact model) is used to consider fluid-structure interaction, combining the pore network model with the particle discrete element model. The main body of the hydraulic fracture propagation particle discrete element numerical algorithm consists of a model initialization module, a flow time-step update module, and a solid time-step update module. The algorithm uses sequential coupling to achieve iterative updates of fluid exchange and particle motion / stress state within the pore network, based on the updated pore domain volume V. d and flow channel opening w p Determine whether the particle connections are broken. If so, the cracks propagate. This process is used to simulate the hydraulic crack propagation process in a homogeneous medium, and numerical simulation results are obtained.

[0024] Discrete element models of conglomerate particles with different gravel characteristics were constructed based on the Voronoi polygon-based gravel generation algorithm. ① Initial Gravel-Matrix Particle Generation: Based on the actual gravel size distribution and volume content, a close-packed aggregate of initial gravel and matrix particles is generated. Disk-shaped particles with the same volume as the gravel are used to approximate the gravel size, characterizing the statistical distribution of gravel size.

[0025] ② Gravel-matrix particle contact topology identification: Based on the polygonal loop algorithm, several particles are connected end-to-end through contact to form a closed loop, constituting a basic topological unit. This closed loop is a convex polygon that does not contain any closely adjacent particles.

[0026] ③ Gravel-Matrix Polygon Particle Topology Identification: Connect the centroids of adjacent closed loops to establish the gravel-matrix polygon particle topology. Each polygon particle is a convex polygon that does not contain any closely adjacent particles, and each convex polygon uniquely corresponds to an initially generated gravel particle or matrix particle. The polygon corresponding to the gravel particle is used as the gravel mesh, and the polygon corresponding to the matrix particle is used as the matrix mesh. Simultaneously, the line-to-line contact between adjacent polygon particles characterizes the geometric features of the linear contact between the gravel-matrix interface and the gravel-gravel interface, and allows for the targeted allocation of contact models and contact parameters between matrix, gravel, and interface particles using polygon boundaries.

[0027] ④ Gravel-Matrix Mesh Refilling: Based on the polygonal particle topology, small-sized disk particles are used to refill the model's computational domain. Disk particles within the gravel mesh region are considered gravel constituent particles, and those within the matrix mesh region are considered matrix constituent particles. In this case, particle contact types can be categorized as matrix-matrix particle contact within the matrix, gravel-gravel particle contact within the gravel, gravel-matrix interface particle contact, and gravel-gravel interface particle contact. Gravel-matrix interface particle contact and gravel-gravel interface particle contact at the gravel-matrix mesh boundary are assigned smooth joint contact models and given corresponding joint surface contact parameters; matrix-matrix particle contact and gravel-gravel particle contact within the matrix and gravel meshes are assigned linear parallel bonded contact models and given corresponding bond parameters.

[0028] Next, the flow time step is updated. S2, Flowing Time Step Update S2-1, Contact Force Update: Based on the particle contact model, the motion state and force state of the rigid particle system are explicitly iteratively updated. In S2-1, the specific method for explicitly iteratively updating the motion state and force state of the rigid particle system is to establish the contact force update equation. At the same time, this step also updates the resultant force of the surrounding pore fluid pressure on a single particle. The gravel-matrix particle contact model has been updated to: , ,

[0029] In the formula, F l and The contact forces, in N, are the linear and parallel bonding portions, respectively. The contact torque of the parallel bonding portion is N·m; It is the unit normal vector perpendicular to the plane of the paper and pointing outwards; and N represents the normal and tangential components of the linear partial contact force. and For the normal and tangential components of the contact force in the parallel bonded area, N; M b The bending moment is N·m. For planar normal contact force, N, The tangential contact force is N; The gravel interface contact model has been updated to...

[0030] In the formula, and These are the normal and tangential contact forces on the joint surface, respectively, in N. Let N be the normal contact force in the plane of any element j. Let N be the tangential contact force of any element j.

[0031] S2-2, Flow channel opening update: Based on the normal contact force on the associated particles, i.e., the force state, update the iterative equation and update the flow channel opening w. p ; In step S2-2, when the contact adhesion is intact and the object is in a compressed state... ; In the formula, w 0 represents the initial opening of the flow channel, in meters (m). F n0 The flow channel opening is w The normal contact force at 0 / 2, in N; when the contact adhesion is intact and under tension, ; When the contact bond breaks ; In the formula, w b The flow channel opening, in meters, corresponds to the instantaneous breakage of the contact bond. m Δ is a dimensionless correction factor for the flow channel opening. w d The change in particle contact gap within the current time step is represented by m, w0, and m, which are variables related to the permeability properties of the pore network, and their values ​​are determined through parameter correction. It is the normal contact force.

[0032] Laminar flow equations based on the pressure difference between the two ends of an incompressible fluid in the pore domain:

[0033] In the formula, q p The fluid volumetric flow rate, m 3 / s; w p Let m be the flow channel opening. p 2 and p 1 represents the fluid pressure in pore domain 1 and pore domain 2, respectively, in Pa; μ The viscosity of the fluid is Pa·s; l p Let be the length of the flow channel, in meters (m). w 0 and m All of these are variables related to the permeability properties of pore networks, and their values ​​are determined through parameter correction.

[0034] S2-3. Obtain the fluid volume change Δ within a single pore domain during a single flow step. V fThe fluid volume change Δ is obtained by measuring the injection amount from the source term and the fluid injection amount from the flow channel. V f ; The change in fluid volume Δ within a single pore domain during a single flow step V f It consists of the source term injection volume and the flow channel fluid injection volume.

[0035] In S4,

[0036] In the formula, The first pore domain connected to this pore domain i The volumetric flow rate of the injected fluid in each flow channel is m³. 3 / s; The first one connected to this pore domain i Fluid volumetric flow rate in each flow channel, m 3 / s; n This represents the total number of flow channels connected to this pore region; Δ t f1 Let be the flow time step, s. The change in pore pressure is influenced by the combined effects of changes in fluid volume and pore volume.

[0037] S2-4. Calculate the pore pressure change Δp in each pore of the pore network; In S2-4, the pore pressure change Δp of each pore in the pore network is specifically as follows: Changes in fluid volume within pores lead to changes in pore pressure, which in turn alters the net flow pressure on the particles, affecting their motion. This change in particle motion, in turn, acts on the pore volume, thus achieving bidirectional coupling. The above process is represented as follows:

[0038] In the formula, Δ p The change in pore pressure is expressed in Pa. K f Here, Δ is the fluid bulk modulus, Pa; V p The change in pore volume is m. 2 ; V p For true pore volume, m 2 .Pick V p = V d φ r ( φ r (Porosity).

[0039] S2-5. Determine the steady flow time step and obtain the critical flow time step Δ.t fc ,

[0040] In the formula, m represents the opening of the i-th flow channel connected to the mother pore. Let m be the length of the i-th flow channel connected to the mother pore. K f The bulk modulus of the fluid, Pa. V p For true pore volume, m 2 , This represents the change in pore volume. φ r Porosity is the fluid viscosity, mPa·s, and n is the number of particles; Introducing a safety factor α t By replacing the first term on the right-hand side of the above equation, we obtain the steady-state flow time step.

[0041] In the formula, 0 < α t <1, take α t It is 0.8; Discrete element method (DEM) updates variables in an explicit format; therefore, each iteration step requires constraint calculations to ensure the numerical stability of the model. Consider having... n The main pore of each flow channel, under the driving force of pressure difference, will cause the fluid inside the pore to interact with the surrounding fluid. n Fluid exchange occurs within the individual pores. Within a single flow step, the direction of the pressure gradient of natural seepage does not change due to fluid exchange; that is, the change in pore pressure in the parent pore cannot exceed the minimum driving pressure difference of the flow channel.

[0042] In the formula, Let Pa be the driving pressure difference of the i-th flow channel connected to the mother pore. Assume... n Each sub-pore injects fluid into the parent pore, while considering Equal and Critical situation The critical flow time step Δ can be obtained. t fc :

[0043] In the formula, m represents the opening of the i-th flow channel connected to the mother pore. Let m be the length of the i-th flow channel connected to the mother pore. Note that the first term on the right-hand side of the above equation is always greater than or equal to 1. A safety factor is introduced. α t By replacing the first term on the right-hand side of the above equation, we obtain the steady-state flow time step.

[0044] In the formula, 0 < α t <1, take α t It is 0.8.

[0045] S3, compared to the steady flow time step Δ t f With stable solid time step Δ t m The smaller of the values ​​is taken as the final iteration step; in actual numerical calculations, the steady-state flow step Δ is compared with the time step Δ. t f With stable solid time step Δ t m The smaller of the two values ​​is used as the final iteration step.

[0046] S4. Solid time step update: Based on the final iteration time step, repeat step S2 for alternating iterative updates. This involves determining the parameters from the previous solid time step, including particle stress state, particle displacement, particle position, and particle contact force. Finally, update the particle velocity and position, and pore volume V of the current time step. d and flow channel opening w p We determine whether the updated particle state meets the particle equilibrium condition. If it does, it means that the crack propagation is complete. This is used to simulate the hydraulic crack propagation process in a homogeneous medium and obtain numerical simulation results. The final iteration time step is determined based on the flow time step module parameters (fluid flow time step Δtf, fluid volume exchange ΔVf, pore pressure change Δp in each pore of the pore network, and the resultant force of the surrounding pore fluid pressure on a single particle). This determines the parameters of the previous solid time step module (particle stress state, particle displacement, particle position, and particle contact force). The particle velocity and position, pore volume Vd, and flow channel opening wp of the current time step are updated. It is then determined whether the updated particle state satisfies the particle equilibrium condition. If so, crack propagation is complete.

[0047] Using the Particle Discrete Element Method (PFC2D) software, Fish program code was written for a conglomerate sample from the Mahu block. A microscale hydraulic fracture propagation model of the conglomerate based on the Voronoi polygon grain-based gravel generation algorithm and fluid flow pore network was established. The dynamic evolution process of hydraulic fractures in the conglomerate was quantitatively analyzed. A systematic numerical simulation study was conducted on the influence of gravel development characteristics in sandstone and conglomerate on the hydraulic fracture propagation law, revealing the interaction mechanism between hydraulic fractures and gravel in conglomerate reservoirs.

[0048] Reference parameters for simulation of hydraulic fracture propagation in conglomerate

[0049] (1) Simulated crack propagation morphology when gravel content is 50% and 90%, results attached. Figure 5 As shown, the gravel content affects the size of the space for free propagation of fractures within the matrix. The lower the gravel content, the less interaction occurs between the fractures and the gravel. For argillaceous conglomerate, when... V g When the crack propagation rate is 50%, the crack propagation path consists of straight cracks within the matrix and tortuous gravel edge cracks. V g When the density is 90%, the gravel is closely adjacent in a linear contact manner, and the fracture lacks free propagation space. However, under the dominant stress, the fracturing fracture seeks to extend along the gravel path with the least energy consumption, and the overall propagation direction is approximately along... σ H For sandy conglomerate, when V g When the crack coverage is 50%, the crack propagation path consists of straight cracks within the matrix, localized cracks around the gravel, and localized cracks penetrating the gravel. V g When the crack morphology is 90%, the post-compression crack morphology is more complex.

[0050] (2) The crack propagation morphology was simulated when the gravel size was 2-8 mm (small conglomerate) and 16-32 mm (large conglomerate). The results are shown in the attached figure. Figure 5As shown, small conglomerates tend to form meandering fractures around the gravel after compression, while large conglomerates exhibit various fracture propagation patterns, including those around the gravel, through the gravel, and intra-gravel arrest. The gravel diameter influences the fracture morphology and type in the near-wellbore area; larger gravel diameters are more prone to forming through-gravel fractures, and fracture propagation in the near-wellbore area is more easily suppressed. With increasing gravel diameter, the number of through-gravel fractures and shear fractures increases, and the fracture extension direction gradually becomes dominated by the gravel. As the gravel diameter increases from 2-8 mm to 16-32 mm, the proportion of around-gravel fractures decreases from 63.6% to 50.9%, while the proportion of through-gravel fractures increases significantly from 15.0% to 49.1%. The proportion of tension fractures in the conglomerate decreases slightly from 99.0% to 98.3%, while the proportion of shear fractures increases slightly from 1.0% to 1.8%. Under the same fracturing scale, with increasing gravel diameter, the initiation and propagation of multi-branch fractures become more difficult, branch fractures are more prone to intra-gravel arrest, and main fractures extend through the gravel to the boundary. The fracture extension direction is dominated by the gravel shape. As the gravel diameter increases, the complexity of near-wellbore fractures decreases. The half-fracture length first increases and then decreases. Medium gravel diameter (9-15mm) is conducive to the formation of a complex fracture network that controls near-wellbore fractures and extends far-wellbore fractures.

Claims

1. A method for predicting crack propagation morphology considering gravel characteristics, characterized in that: Includes the following steps, S1. Based on the Voronoi polygonal grain-based gravel generation algorithm, construct discrete element models of conglomerate particles with different gravel characteristics. Assign smooth joint contact models to the gravel-matrix interface and the gravel-gravel interface, and assign corresponding joint surface contact parameters. Assign linear parallel bonded contact models to the matrix-matrix particle contact and the gravel-gravel particle contact, and assign corresponding bonded parameters. S2, Flowing Time Step Update S2-1, Contact Force Update: Based on the particle contact model, the motion state and force state of the rigid particle system are explicitly iteratively updated. S2-2, Flow channel opening update: Based on the normal contact force on the associated particles, i.e., the force state, update the iterative equation and update the flow channel opening w. p ; S2-3. Obtain the fluid volume change Δ within a single pore domain during a single flow step. V f The fluid volume change Δ is obtained by measuring the injection amount from the source term and the fluid injection amount from the flow channel. V f ; S2-4. Calculate the pore pressure change Δp in each pore of the pore network; S2-5. Determine the steady flow time step and obtain the critical flow time step Δ. t fc , In the formula, m represents the opening of the i-th flow channel connected to the mother pore. The length of the i-th flow channel connected to the mother pore is m; K f Here, denoted as the fluid bulk modulus, is expressed in Pa. V p For true pore volume, m 2 , This represents the change in pore volume. φ r Porosity is the fluid viscosity, mPa·s, and n is the number of particles; Introducing a safety factor α t By replacing the first term on the right-hand side of the above equation, we obtain the steady-state flow time step. In the formula, 0 < α t <1, take α t It is 0.8; S3, compared to the steady flow time step Δ t f With stable solid time step Δ t m The smaller of the two values ​​is taken as the final iteration step; S4. Solid time step update: Based on the final iteration time step, repeat step S2 for alternating iterative updates. This involves determining the parameters from the previous solid time step, including particle stress state, particle displacement, particle position, and particle contact force. Finally, update the particle velocity and position, and pore volume V of the current time step. d and flow channel opening w p We determine whether the updated particle state satisfies the particle equilibrium condition. If so, it indicates that the crack propagation is complete. This is used to simulate the hydraulic crack propagation process in a homogeneous medium, and numerical simulation results are obtained.

2. The method for predicting crack propagation morphology considering gravel characteristics according to claim 1, characterized in that: In S1, the specific steps for constructing discrete element models of conglomerate particles with different gravel characteristics are as follows: Gravel-matrix initial particle generation: Based on the actual gravel size distribution and volume content, a close-packed mass of gravel and matrix initial particles is generated; Gravel-matrix particle contact topology identification: Based on the polygonal loop algorithm, several particles are connected sequentially through the contact head and tail to form a closed loop, which constitutes the basic topological unit; Gravel-matrix polygonal particle topology identification: Connect the centroids of adjacent closed loops to establish the gravel-matrix polygonal particle topology; Gravel-matrix mesh region refilling: Based on the polygonal particle topology, small-sized disk particles are used to refill the model's computational domain. Disk particles located within the gravel mesh region are considered as gravel constituent particles, and disk particles located within the matrix mesh region are considered as matrix constituent particles. At this point, particle contact types can be classified into matrix-matrix particle contact within the matrix, gravel-gravel particle contact within the gravel, gravel-matrix interface particle contact, and gravel-gravel interface particle contact. Smooth joint contact models are assigned to the gravel-matrix interface particle contact and gravel-gravel interface particle contact at the gravel-matrix mesh boundary intersection, and corresponding joint surface contact parameters are assigned; linear parallel bonded contact models are assigned to the matrix-matrix particle contact and gravel-gravel particle contact within the matrix mesh and gravel mesh, and corresponding bonded parameters are assigned.

3. The method for predicting crack propagation morphology considering gravel characteristics according to claim 1, characterized in that: In S2-1, the specific method for explicitly iteratively updating the motion and force states of the rigid particle system is as follows: The gravel-matrix particle contact model has been updated to: ， ， In the formula, F l and The contact forces, in N, are the linear and parallel bonding portions, respectively. The contact torque of the parallel bonding portion is N·m; It is the unit normal vector perpendicular to the plane of the paper and pointing outwards; and N represents the normal and tangential components of the linear partial contact force. and For the normal and tangential components of the contact force in the parallel bonded portion, N; M b The bending moment is N·m. For planar normal contact force, N, The tangential contact force is N; The gravel interface contact model has been updated to... In the formula, and These are the normal and tangential contact forces on the joint surface, respectively, in N. Let N be the normal contact force in the plane of any element j. Let N be the tangential contact force of any element j.

4. The method for predicting crack propagation morphology considering gravel characteristics according to claim 1, characterized in that: In step S2-2, when the contact adhesion is intact and the object is in a compressed state... ; In the formula, w 0 represents the initial opening of the flow channel, in meters (m). F n0 The flow channel opening is w The normal contact force at 0 / 2, in N; When the contact bonding is intact and the object is under tension ; When the contact bond breaks ; In the formula, w b The flow channel opening, in meters, corresponds to the instantaneous breakage of the contact bond. m This is a dimensionless correction factor for the flow channel opening. Δ w d Let m be the change in particle contact gap within the current time step. w 0 and m All of these are variables related to the permeability properties of pore networks, and their values ​​are determined through parameter correction. It is the normal contact force.

5. A method for predicting crack propagation morphology considering gravel characteristics according to claim 4, characterized in that: In the formula, q p The fluid volumetric flow rate, m 3 / s; w p Let m be the flow channel opening. p 2 and p 1 represents the fluid pressure in pore domain 1 and pore domain 2, respectively, in Pa; μ The viscosity of the fluid is Pa·s; l p Let be the length of the flow channel, in meters (m).

6. The method for predicting crack propagation morphology considering gravel characteristics according to claim 1, characterized in that: In S4, In the formula, The first pore domain connected to this pore domain i The volumetric flow rate of the injected fluid in each flow channel is m³. 3 / s; The first one connected to this pore domain i Fluid volumetric flow rate in each flow channel, m 3 / s; n This represents the total number of flow channels connected to this pore region; Δ t f1 Let s be the time step during the flow.

7. The method for predicting crack propagation morphology considering gravel characteristics according to claim 1, characterized in that: In S2-4, the pore pressure change Δp of each pore in the pore network is specifically as follows: In the formula, Δ p The change in pore pressure is expressed in Pa. K f Here, Δ is the fluid bulk modulus, Pa; V p The change in pore volume is m. 2 ; V p For true pore volume, m 2 .

Images